co 

■^r 


<£ 

i/n 



LOAN COPY: RET ^ ms 
AFVVL (WLlL^J^g 
KIRTLAND AFB, f wig 


DETAILED DESCRIPTION AND RESULTS 
OF A METHOD FOR COMPUTING MEAN 
AND FLUCTUATING QUANTITIES IN 
TURBULENT BOUNDARY LAYERS 


by Ivan E. Beckwith and Dennis M. Bushnell 

Langley Research Center , < v 

• :* iV’ 1 ' 

Langley Station , Hampton , Va. Jp * \ 

i ■ 

• £ " 

NATIONAL AERONAUTICS AND SPACE ADMINISTRATION • WASHINGTON, D. C. . OCTOBER 1968 


BRARY KAFB, NM 



TECH LIBRARY KAFB, NM 



□131715 


DETAILED DESCRIPTION AND RESULTS OF A METHOD 

FOR COMPUTING MEAN AND FLUCTUATING QUANTITIES IN 

TURBULENT BOUNDARY LAYERS 

By Ivan E. Beckwith and Dennis M. Bushnell 

Langley Research Center 
Langley Station, Hampton, Va. 


NATIONAL AERONAUTICS AND SPACE ADMINISTRATION 


For sale by the Clearinghouse for Federal Scientific and Technical Information 
Springfield, Virginia 22151 — CFSTI price $3.00 




DETAILED DESCRIPTION AND RESULTS OF A METHOD 
FOR COMPUTING MEAN AND FLUCTUATING QUANTITIES IN 
TURBULENT BOUNDARY LAYERS 

By Ivan E. Beckwith and Dennis M. Bushnell 
Langley Research Center 

SUMMARY 

The conservation equations for mass, mean momentum, and turbulent kinetic energy 
for the incompressible turbulent boundary layer have been solved by an implicit finite- 
difference procedure. Mathematical models developed by Glushko (Bull. Acad. Sci. USSR, 
Mech. Ser., no. 4, 1965) for the production, dissipation, and diffusion of the turbulent 
kinetic energy in the flat- plate turbulent boundary layer have been modified and used to 
calculate a nonequilibrium turbulent boundary layer subjected initially to a large adverse 
pressure gradient followed by a run of constant pressure. 

Comparisons of both mean and fluctuating flow properties have indicated generally 
good agreement between the calculated results and experimental measurements of 
Goldberg (MIT Rep. No. 85, 1966). The best overall agreement with data was obtained by 
reducing the scale of turbulence in the outer part of the boundary layer to about 70 percent 
of the flat-plate values as used by Glushko. The calculations have indicated that further 
simple modifications to the turbulence scale function and to some of the mathematical 
models for the turbulence correlation terms should improve the accuracy of predictions 
for the Goldberg data. The present authors have shown (paper presented at 1968 Confer- 
ence on Computation Methods in Turbulent Boundary Layers at Stanford University, 

Aug. 1968) that predictions in good agreement with data were obtained for other arbitrary 
pressure distributions. 


INTRODUCTION 

The basic assumption in conventional methods for the calculation of turbulent flows 
is that the mean values of products of the fluctuating characteristics of the flow are 
related in some unique manner to the mean flow properties. Such relations between 
these correlations of fluctuating flow quantities and the mean flow cannot be determined 
solely from theoretical analysis. In the Reynolds equations of mean motion, this problem 
of relating the correlations to the mean flow has generally been met by introducing more 
or less arbitrary assumptions for the Reynolds stress terms (ref. 1, pp. 277-293). In 
1945, Prandtl, Nevzgljadov, and Chou (refs. 2 to 5) reported on independent investigations 



intended to develop a more rigorous approach to this problem. These investigations were 
based on the idea of using independent differential equations to describe the dynamics of 
the correlations for the turbulent velocity fluctuations. The correlation equations are 
derived (ref. 1, pp. 250-260, for example) from the Navier-Stokes equations of motion and 
contain terms for double velocity correlations (Reynolds stress and turbulence kinetic 
energy terms), triple velocity correlations, and correlations of velocity and pressure 
fluctuations. 

In 1951 Rotta (ref. 6) extended the work of Prandtl (ref. 2), and, in particular, that 
of Chou (refs. 4 and 5), in considerable detail; his work was based on advances in theory 
and new data not available during the earlier investigations of 1945. (See also discussion 
in ref. 7, pp. 43-54.) These methods were developed further and applied to various types 
of simple flows as described in references 8 to 11. An integral form of the turbulence 
energy equation was utilized by McDonald (ref. 12) to compute the mean profiles by an 
integral method. 

Kovasznay and Nee (refs. 13 and 14) assumed that the effective total viscosity obeys 
a "rate equation" expressing the "conservation" of the total viscosity. Harlow and 
Nakayama (ref. 15) gave a more formal derivation (based on the equation for the turbu- 
lence kinetic energy) of a rate equation for the eddy viscosity and also constructed a 
"transport" equation for the scale of turbulence by analogy with Brownian motion. The 
same authors (ref. 16) have since derived a new transport equation for the dissipation 
function. Since this function depends on the scale of turbulence, their new equation and 
their previous rate equation for the eddy viscosity were proposed as the basic equations 
for a general method of computing turbulent shear flows. 

In all the investigations mentioned, except the integral method of McDonald, assump 
tions regarding the relative magnitudes of the various fluctuating quantities were made in 
order to simplify and obtain solutions to the nonlinear equations involved. Consequently, 
these methods were applied mostly to simple flows or portions of simple flows where 
some terms in the correlation equations could be neglected. Hence, the mathematical 
formulations of remaining terms representing the fluctuating quantities could not be 
tested for their degree of generality. In the integral method of McDonald, it is still not 
possible to test the detailed spatial variations of these formulations and the turbulent dif- 
fusion terms drop out completely upon integration across a shear layer. 

The answer to these difficulties is, of course, to obtain numerical solutions of the 
complete equations with automatic computing machines. Such solutions have been given 
recently by Glushko (ref. 17), Bradshaw, Ferriss, and Atwell (ref. 18), and Nash (ref. 19). 
In the methods of references 18 and 19, the molecular shear is neglected; it was therefore 
necessary to provide the correct wall boundary condition by incorporating the "law-of-the 
wall" relation between velocity and wall shear. Glushko, on the other hand, kept all the 
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viscous terms and used formulations with essentially correct limiting forms at the wall. 

It might be expected that this latter approach would therefore be somewhat more general 
than that of Bradshaw in that sudden changes in wall-boundary conditions could be negoti- 
ated and the extension to compressible flows where the law-of-the-wall relation may not 
be generally applicable should give better results. The only computation presented by 
Glushko, however, was for the flat plate; again, the degree of generality of his assumptions 
for the fluctuating flow parameters could not be determined. 

The primary purpose of the present paper is to test the method of Glushko in a non- 
equilibrium adverse-pressure- gradient flow and to determine whether his formulations of 
the turbulence quantities based primarily on flat-plate data result in satisfactory predic- 
tions for the mean properties of this more complex flow. Since the ultimate success of 
these methods depends on the assumptions for the fluctuating properties, comparisons of 
values from numerical solutions for these quantities with experimental data are also 
made. Except for shear profiles in references 14, 18, and 19, and turbulent kinetic energy 
profiles in reference 17, such comparisons for boundary-layer flows have not been pub- 
lished yet. The work reported herein is part of a general investigation intended to develop 
these methods for application to compressible flows. 

This paper includes an appendix by Carolyn C. Thomas, of the Langley Research 
Center. This appendix presents a description and listing of the digital computer program. 

SYMBOLS 

A n ,B n ,C n ,D n functions in linearized finite-difference form of momentum equation 
(eq. (B13)) 

A n ,B n ,C n ,D n functions in linearized finite-difference form of energy equation 
(eq. (B27)) 

C constant in model of turbulent dissipation term (eq. (12)) 


C 1>C2, C 3,C4, C 5,C^1 

C’,Ci’,C2',C 3 ’,C 6 ’ J 


constants used in appendix A 


Cf 


local skin friction, 
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D 


dimensionless diffusion and dissipation function 1 + Q!kRH(kR) (eq. (19)) 


E turbulent kinetic -energy profile function, e /u e 2 


E , 


constant in expression used for initial E profile (eq. (31a)) 


1 / 2 2 2 \ 

instantaneous turbulent kinetic energy per unit mass, — (u^' + u 2 * + u 3 * / 


1/2 2 2 

mean turbulent kinetic energy per unit mass, —I uj T + ug' + ug' 


F 

f 

GruSn 

H(R) 

H* 


K 


K = 2 


1 + K 


_ 2K 


1 + K 


l 

M 

m,n 

N 

n 


P 

4 


velocity profile function, uiyUe 

factor in expression used for flat-plate solutions (see expression (31b)) 

functions required to calculate F m+ i (eq. (B16)) 
eddy viscosity function (eq. (9b)) 
shape factor, 5*/ 6 

factor for variable A 77 step size, A 7 ] n /Arj n _i 


exponent in A£ expression used for flat -plate solutions 
(see expression (31b)) 

characteristic mean scales of turbulent motion 
dimensionless total viscosity function 1 + oRH(R) (eq. (19)) 

computing grid location (see sketch in appendix B) 

maximum value of n 

exponent in definition of 77 (eq. ( 22 )) 

pressure 



fluctuation in pressure, p - p 
mean pressure 


yet 

turbulent Reynolds number, 

constant in models of turbulent terms (eq. (9)) 

u e x l 

Reynolds number based on xj, — - — 

U e 6 

Reynolds number based on 6, — — 

U e 0 

Reynolds number based on 6, — — 
time 

local free-stream velocity in xj direction at edge of boundary layer 
instantaneous velocity in x* direction 
mean velocity in x* direction 
fluctuation in velocity in x^ direction, Uj - Uj 
friction velocity, ^ T w /P 
transformed normal velocity (eq. (26)) 

Cartesian coordinates in tensor notation (i = 1,2,3) 
constant in models of turbulent terms (eqs. (9) and (10)) 
boundary-layer thickness ^taken at some specified value of 



displacement thickness, 



dx 2 


eddy viscosity (eq. (7)) 



e e error criteria, 1 - F e 

e e ' error criteria (eq. (B18)) 

e w convergence criteria (eq. (B24)) 

f 6 /u i uj 2 \ 

6 momentum thickness \ dx 2 

J 0 \ u e u ej 

k constant in models of turbulent terms (eqs. (12) and (13)) 

ju dynamic viscosity 

v kinematic viscosity, \i/p 

£,77 similarity coordinates (eqs. (21) and (22)) 

77* constant in expression used for initial E profile (eq. (31a)) 

p density 

9ui 

r local shear stress, p puj'u2' 

9x 2 

Trj, turbulent shear stress, -puj'^' 

4> turbulent scale function (eq. (8)) 

Subscripts: 

average value 
dissipation 

Ul 

edge of boundary layer, = (*■ “ e e) 

U e 

standard tensor notation: 1, 2, and 3 denote directions that are, respec- 
tively, parallel to the surface and in the same direction as the external 
velocity vector, normal to the surface, and parallel to the surface but 
normal to the external velocity vector 


av 
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o initial conditions except where noted 

T turbulent 

t transition 

w wall or surface 

6 evaluated at some specified value of F 

REVIEW OF PROBLEM 

The following remarks are not intended as a complete or general review of the 
literature on turbulent boundary -layer theory. Rather, it is desired to emphasize certain 
aspects of previous work that are of particular significance for the present extension of 
the Glushko method (ref. 17). 

Applications of Turbulent Kinetic Energy Equation 

From physical and dimensional considerations, Prandtl (ref. 2) derived the equation 
for the conservation of kinetic energy of turbulent velocity fluctuations e (one of the 
double velocity correlations) and compared the calculated variation of e across a chan- 
nel with experimental data. Nevzgljadov (ref. 3) presented a more elegant derivation of 
the equation for e and applied it, together with the equations for the mean motion, to the 
flow in a circular pipe. In these methods it was necessary to assume functional relations 
between e and the various other fluctuation terms which appear, such as the production 
and dissipation of turbulent energy; an "eddy" viscosity concept was used to formulate the 
Reynolds stress. (The dominant production term is usually the product of the Reynolds 
stress and the normal gradient of the mean velocity.) The success of the methods depends 
on how well the assumed functional relations represent the real situation. A character- 
istic length l which is interpreted as the mean scale of the turbulence enters some of the 
functions, In reference 3 the assumption that l °c e was used; whereas in reference 2 it 
was assumed that l was related directly to the mixing length and, hence, depended on 
the type and geometry of the flow. In spite of these different assumptions for l, the 
agreement with experiment shown in both references 2 and 3 was reasonable. This agree- 
ment indicates that assumptions relating the different fluctuating properties (actually, 
mean values of products of fluctuating quantities) of the flow to each other may not be as 
critical as assumptions relating fluctuating flow properties directly to mean flow 
quantities. 
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In the method proposed by Chou (ref. 4) and applied by him to channel and flat-plate 
flows (refs. 5 and 20), separate equations for each of the different components of both the 
double and triple velocity correlations were used and the concept of an eddy viscosity was 
not required. If quadruple correlations are neglected, this set of equations then repre- 
sents a closed system that could be used to compute any turbulent flow if sufficiently gen- 
eral formulations for terms that represent the diffusion and dissipation of turbulence can 
be determined. 

Emmons (ref. 8) further developed and applied the method of Prandtl (ref. 2) to the 
calculation of various turbulence quantities for channel, free jet, and wake flows. By 
appropriate adjustment of the constants and the turbulence scale factor l, good agree- 
ment with corresponding measurements (refs. 21 to 23) was obtained. Townsend (ref. 9) 
assumed that the turbulent shear is directly proportional to e (and thereby abandoned 
the eddy viscosity concept) within the inner region of a shear layer and used the turbulent 
energy equation to compute mean velocity distributions for equilibrium shear flows. 

Levin (ref. 10) applied Rotta’s method to pipe flow and obtained reasonable agreement 
with the experimental data of Laufer (ref. 24) for turbulent shear and the production and 
dissipation of turbulent energy. Spalding (ref. 11) has applied the turbulence energy 
equation to separated flows and used assumptions for the turbulence quantities similar to 
those of Prandtl (ref. 2). 

As mentioned in the Introduction, numerical solutions to the complete equations are 
essential to determine the degree of generality of the models used for the fluctuating 
terms. Although numerical procedures were used in references 14 and 15, solutions 
were given only for simplified forms of the equations applicable to flat plate or pipe flows. 
In the solutions of references 17 to 19, the complete equations were retained except for 
the neglect of molecular shear in references 18 and 19. The equations of mean motion, 
continuity, and turbulence energy then become hyperbolic (ref. 18) and can therefore be 
solved by the method of characteristics. The numerical computing procedures are 
thereby considerably simplified but the wall boundary condition of zero turbulent shear 
cannot be satisfied directly. Nevertheless, Bradshaw, Ferriss, and Atwell (ref. 18) 
obtained good agreement of mean velocity, shear profiles, and the variation of thickness 
parameters with experimental data for both equilibrium and some nonequilibrium adverse 
pressure gradient cases. Bradshaw and Ferriss have applied the method of reference 18 
to compressible flows (ref. 25) at low supersonic Mach numbers with only partial success. 
The method of Nash (ref. 19) differs from that of reference 18 only in that the computing 
procedure is a more conventional finite-difference method rather than the method of 
characteristics as applied uniquely by Bradshaw. 
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Integral Methods 


In the conventional integral methods, the mean momentum equation is integrated 
across the boundary layer and solutions are then obtained by the use of assumptions for 
velocity profiles (or auxiliary equations for the shape factor, H* = 6 */ 8 ) and skin friction 
as functions of 6 or other mean flow parameters. Large discrepancies in predicted 
variations of H*, 6 , and Cf result from these methods, as shown by the recent review 

of Thompson (ref. 26). In another integral method that has been used extensively, the 
integral form of the mean kinetic energy equation is solved simultaneously with the mean 
flow integral momentum equation as originally done by Truckenbrodt (ref. 27). The addi- 
tional requirement here is to develop a sufficiently general relation for the dissipation 
integral. A recent correlation scheme for this integral as described by Walz (ref. 28) 
has resulted in good agreement with experimental distributions of H* in several 
adverse-pressure-gradient cases. The dissipation integral depends on the turbulent 
shear distribution across the boundary layer; however, the various integral methods pro- 
vide no information on local shear or other local turbulence parameters and hence can be 
evaluated only in terms of mean flow properties such as H*, 9 , and Cf. Although the 

auxiliary relations for H*, Cf, and the dissipation integral can be tailored to give agree- 
ment with a specific set of data, extension to other situations, particularly compressible 
flow, always remains doubtful. 

Finite-Difference Procedures 

An approach which is more basic than the integral methods is to solve the mean 
momentum and continuity equations directly by finite-difference procedures so that the 
only unknown function is the relation between the Reynolds stress and mean flow proper- 
ties. The success of this approach as applied to several different types of boundary 
layers (refs. 29 and 30, for example) indicates that the difficulties in the integral methods 
(besides the basic limitations of any boundary -layer integral method) are caused mainly 
by deficiencies in the assumptions for the velocity profiles (or H*), the shear stress and 
dissipation integrals, and the skin-friction relations. In the finite -difference methods of 
references 29 and 30, simple functions for the eddy viscosity e in terms of local flow 
conditions are used. The Reynolds stress is thus assumed to depend directly on the local 
mean flow velocity gradient and other local mean flow properties. 

On the other hand, in the finite -difference procedures of references 17 to 19, the 
Reynolds stress is also related to correlations of fluctuation quantities by the turbulent 
kinetic energy equation. Formulations for the turbulence correlation quantities in ref- 
erence 17 were based primarily on the work of Prandtl and Rotta (refs. 2 and 6) where 
an eddy viscosity concept is used. Assumptions for the correlations in reference 18 are 
based on extensions of the ideas of Townsend (ref. 9); that is, all turbulence quantities 
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were related directly to the turbulent shear stress and the turbulence energy equation 
becomes an equation for the shear stress. These methods combine the advantages of 
a finite -difference solution of the mean momentum equation with the more general 
relation for the turbulent shear as supplied by the turbulence kinetic energy equation. 

BOUNDARY- LAYER EQUATIONS 


For conditions where a characteristic Reynolds number \5 e 6/v is at most of the 
order of xj^6 and if x-^j 6 » 1, the two-dimensional time-steady equations of motion 
for the turbulent flow of an incompressible fluid (ref. 1, pp. 453-457) reduce to: 

For xj mean momentum: 


Ul 


0U] 

0Xi 


+ u 2 


9uj 

9Xr 


l^e 

P dx n 


9 

9Xr 


Uj u 2 + v 


S^Uj 

9Xo' ! 


( 1 ) 


For turbulent kinetic energy: 


- 9e 

Ul 9x7 


- 9§ 

+ U 2^ 


= -Ui u 


1 u 2 


— 9£L 

t 1 

9xo 


9xc 


u 2 




+ v ■ 


Bh_ 

9x 0 ^ 


9 Ui ’ 0UJ* 


9x i 


8x i 


( 2 ) 


(In accordance with conventional notation, when a subscript i or j is repeated, sum- 
mation over three directions is implied.) These equations together with the equation of 
continuity 


9uj 

9x^ 


9fl 9 

— = 0 
9x 2 


( 3 ) 


and appropriate initial and boundary conditions determine the three dependent variables 
iij, U 2 > and e within the turbulent boundary layer. 

Since the Reynolds shear stress 

T X r ~P u l u 2 (4) 

the diffusion of total turbulence energy by velocity fluctuations 


Diffusion = 


0x 2 


Ur 


— + e 


(5) 
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and the dissipation of turbulence energy 


Bui' W 

Dissipation = v — — (6) 

8 Xj 8 Xj 

are in general unknown, these equations do not represent a closed system. Nevertheless, 
the dependence of the Reynolds shear stress on both the correlation terms and the mean 
flow properties themselves is clearly evident from equations (1) to (3). 

ASSUMPTIONS FOR FLUCTUATING QUANTITIES 

In order to formulate realistic expressions for the correlation of fluctuation terms 
in equation (2) involving Reynolds stress, diffusion of total turbulent energy, and dissipa- 
tion of turbulent energy (eqs. (4) to (6)), it is desirable to rely on experimental measure- 
ments of these quantities whenever possible. Except for the diffusion of pressure energy, 
these terms have been measured in boundary layers with zero pressure gradient 
(Townsend, Klebanoff, and Corrsin and Kistler, refs. 31 to 33) and with pressure gradients 
for both equilibrium (Bradshaw, ref. 34) and nonequilibrium boundary layers (for example, 
Schubauer and Klebanoff (ref. 35) and Bradshaw and Ferriss (ref. 36)). 

The expressions for correlations as developed by Glushko (ref. 17) were based on 
the general approach of Rotta (ref. 6) wherein the dissipation and diffusion terms are 

— p7 

assumed to be functions of e, l, and R = The form of these functions depends pri- 
marily on physical and dimensional reasoning. The Reynolds stress was related to the 
mean velocity gradient by 


Tip — e 


3uj 

9X2 


(7) 


where the dimensionless eddy viscosity e/ju is assumed to be a function only of R 
rather than of the local mean flow properties as in more conventional approaches, such 
as the methods of references 29 and 30. The mean scale of turbulence l was evaluated 
from flat-plate data for two-point correlation coefficients of the longitudinal velocity 
fluctuations in such a way that l is analogous to an integral scale (see ref. 1, p. 37) 
parameter. The corresponding ratio Z/6 was taken as a universal function of the form 
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As indicated by Glushko, there are few data available for the evaluation of l, even for 
flat-plate flows; thus, the generality of this relation remains questionable. However, in 
a recent experimental investigation (ref. 37) of a separating boundary layer, measure- 
ments of the two-point correlation coefficients were obtained. If these measurements 
are used to compute l as defined by Glushko, the resulting 1/5 function is nearly 

x 2 

identical to that of Glushko's for -f- < 0.6. 

6 

With l thus determined, Glushko (ref. 17) then evaluated the function e/p from 
flat-plate data of references 31 and 32 for R > 400 and from the wake data of Townsend 
(ref. 23) for R < 100. Although Glushko was able to obtain reasonable agreement with 
experimental data for e, mean velocity profiles, and Cf on flat plates by adjusting con- 
stants in these formulations, no comparisons with other experimental measurements were 
published. Also, no calculations were reported for adverse-pressure-gradient cases 
which generally provide a more critical test of any theoretical method than the simpler 
flat-plate flows. Before proceeding to compare theoretical predictions with experimental 
data for both mean (including rate of boundary -layer growth) and fluctuating quantities for 
flat-plate and nonequilibrium, adverse -pressure-gradient flows, Glushko’s formulations 
are considered in more detail. 


Reynolds Stress 


The starting point of Glushko’s (ref. 17) derivation of expressions for the fluctuation 
quantities was a relation between the eddy viscosity e/p and the turbulent Reynolds 
number R. This relation had been pointed out previously by Levin (ref. 10) and was 
based primarily on the semiempirical results of Rotta (ref. 6) that, in turn, depend partly 
on the observed decay rates of homogeneous turbulence (ref. 38, pp. 101-113). From flat- 
plate and wake data with l determined from equation (8), as mentioned above, Glushko 
then obtained the simple function 


where 


± = H(R)oR 


(9a) 




<*- < 
Rn 


< JL< 

Rq 




< 


R_<~\ 

R 0 J 


(9b) 
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and a and R 0 are constants. Since this result is based partly on the assumption that 
the dissipation rate of turbulent energy in a boundary layer is similar to that in isotropic 
homogeneous turbulence, it cannot be regarded as firmly established, particularly in the 
very near wall region (where decay of e is caused by the presence of the wall) nor in the 
outer portions of the boundary layer (where the decay in e is associated with the inter- 
mittent character of the flow and the edge boundary condition on e). Indeed, one justifi- 

9 9S 1 

cation for the eddy viscosity concept itself, generally expressed as e = , is perhaps 

8x 2 

that the concept has generally given useful results. However, as shown by Townsend 
(ref. 9) and by Bradshaw, Ferriss, and Atwell (ref. 18), reasonable results can be obtained 
from simultaneous solutions of equations (1) and (2) without the use of the eddy viscosity 
concept. 

Glushko’s final expression for the Reynolds stress is then 


r T = 


-puj’u^ = juH(R)qR 


8uj 

8x 2 


( 10 ) 


For boundary-layer flows, the production of turbulent energy, which is represented by the 
first term on the right in equation (2), is generally large in magnitude, and therefore the 
formulation of t,j, is of critical importance. 


Dissipation 

It is known that the dissipation of energy in homogeneous turbulence is proportional 

to (e) 3//2 /z (ref. 38, p. 106) which has been used directly by Prandtl (ref. 2), Emmons 
(ref. 8), and Townsend (ref. 9) to compute the production of turbulent energy from equa- 
tion (2) for simple shear flows. Based on Rotta's work, Glushko showed that the dissipa- 
tion term in equation (2) can be written as 


v 


9uj' 9m T s 

— - = vC — + 

9 Xj 9 Xj j2 


kC 


(e) 3 / 2 


= vC(l + kR)4 
l 2 


( 11 ) 


and it is seen that the second term on the right is the same as the dissipation in homoge- 
neous turbulence. The first term on the right is important only for small R and whether 
this formulation is applicable to regions of the boundary layer where R is small cannot 
be determined except by comparison with data. Since for large R, the expression already 
obtained for e/p (eqs. (9)) had the desired dependence on R, Glushko assumed the final 
relation for the dissipation term as 


13 



( 12 ) 


8u^' 8u^' 



where — (kR) indicates the same function as given by equations (9) except that R is 
M 

replaced by kR. The turbulence scale l as used by Glushko was considered as the 
same mean scale function given by equation (8); however, as written in equation (12), l 
could be considered a microscale (see ref. 1, p. 37) if C were adjusted accordingly. 
The possibility of introducing a microscale for l in the dissipation term is considered 
in the section "Results and Discussion." 


The dissipation term is generally of the same order of magnitude as the production 
term, except very near the wall, and hence is also of critical importance. At the wall 
itself, equation (12) does not give the correct limiting form for the dissipation which, from 
equation (2), should be 


v 


8uj' 9uj' 
8x. dXz 



An analysis similar to that of Rotta (ref. 7, p. 59) shows that 


{ dx 2 2 


is generally finite, 


w 


whereas equation (12) evaluated at the wall is zero since e w = 0. The use of equation (12) 


thus forces 


8 2 e 

t 8xo 2 


to be zero, which is physically incorrect and has apparently caused 


w 


some numerical difficulties in the present solutions. These difficulties have been largely 
avoided by rejecting negative values of e which sometimes occurred during iteration 
cycles in the numerical procedure. 


Diffusion 

Glushko reasoned that the total diffusion of turbulence energy was due to the gradient 
of e. In addition, he specified the corresponding diffusion coefficient to be the same 
quantity given in brackets in equation (12). Hence, he obtained for the diffusion terms 


8 

8 x 2 
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~ U 2 'f— + el 
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8e { 
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Again, it is obvious that the generality of this assumed expression for the diffusion can 
only be determined by comparison with data. Possible deficiencies are that the diffusion 
of energy by pressure fluctuations is not treated separately and that a simple gradient- 
type diffusion as postulated does not always occur (ref. 1, pp. 288-289). Nevertheless, 
for many practical cases the diffusion terms are smaller than the dissipation and produc- 
tion terms; thus, for these situations, the formulation of diffusion is presumably not 
critical. 


TRANSFORMATION TO SIMILARITY COORDINATES 


Profile Functions 


The main reason for transforming to similarity type of coordinates is to provide 
scale factors that in terms of the transformed variables £ and rj, reduce or remove 
the rate of increase in boundary-layer thickness with distance £ along the surface. The 
number of computing steps A77 required across the boundary layer to obtain desired 
accuracy in the finite-difference procedures can thereby be reduced and kept more nearly 
constant. Also, in a region of approximate local similarity, the streamwise step size A£ 
can be increased since for this situation the rate of change of the dependent variables with 
£ is much reduced. 


Before the boundary -layer equations (eqs. (1) to (3)) are transformed to the £,17 
coordinates, it is convenient to write them in terms of the dimensionless velocity and 
energy profile functions 





(14) 


Then with the use of equations (8) to (13) and Bernoulli's equation for the pressure 
gradient 


U e 


dUe 

dxj 


1 d£e 

P d x 1 


(15) 


equations (1), (2), and (3) become, respectively, 
for mean momentum: 

9F ^2 8F _ 1 - F 2 dU e v 8 f 9f\ 

0X 1 U e 8 x 2 U e dxj U e 9x2 1 8x 2 / 


(16) 


15 



for turbulent kinetic energy: 


F 



u 2 8E 
+ U e 9x 2 


FE du e 
U e dxj 




C 


v DE 



(17) 


and for continuity: 


8F 8 ^2 F dU e 
8xj 9x 2 U e Ue d*! 


(18) 


where 


M = 1 + oRH(R) = 1 + aH(R) <p\[ER 5 
D = 1 + <3«RH(kR) = 1 + Q!kH(/cR) (p\j ERg 


(19) 


The function H(/cR) is identical to the function H(R) (eq. (9b)) except that R/R 0 is 
replaced by kR/R 0 . The first three terms in equation (17) account for the convection of 
turbulent energy and the last three terms represent, respectively, the production, diffusion, 
and dissipation of the turbulent energy. 

The boundary conditions to be applied to this system of equations are 


for x 2 = 0: 


for x 2 — °°: 


or 



F - 1.0 
E - 0 > 

E — E e 


(20a) 


(20b) 


where E e is one-half the square of the free-stream turbulence intensity. 


Similarity Variables 

The transformed variables are defined as 



(21) 
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77( X 1,X 

v U e 

2) = - x 2 
Z7(2|) n 

(22) 

where ii is, in general, a variable function of | and is determined from the require- 
ment that 7?g (the boundary-layer thickness in the transformed coordinates) is constant. 

Hence, from equation (22) 



l°f?e 

n — 

L u e,o u o J 

(23) 

11 — 

log e H 


In a boundary layer where the velocity profile shapes and thickness 6 are changing 
rapidly, the value of n would also change rapidly so that a lengthy iteration procedure 
would be required to satisfy equation (23). An alternate procedure is to specify n as a 
constant; then, for U e constant = UeXj/i'j and 77 6 constant, n = -i for laminar 

boundary layers, and n ~ 0.8 to 1.0 for turbulent boundary layers. (Ref. 39 (p. 537) 
gives h = 0.8 for turbulent boundary layers.) 


Transformed Equations 

The general transformation formulas from equations (21) and (22) may be written as 


8 


0 X 


1 


8 

0x 2 


Ue 8 877 8 

v 8£ 8xj 877 

> 


U. 


v(20 


n 877 


(24) 


The operator that appears on the left-hand side of equations (16) and (17) then becomes 





FU e _8_ | Ue v _8_ 
v 9 £ + K2|) 2S 917 


(25) 


where V is a transformed normal velocity defined as 

V = ( 2 ^) 2 ^ F -?^L + (2|)^ ^ 

U e U e 


(26) 
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Equations (16), (17), and (18) can then be written for momentum: 


(2|) 2fi F n. + V M ^ (l - F 2 ) 


9£ Sr] U p d£ 


Sr] \ 9?7 


(27) 


for energy: 


(2|)25f « + y « . . ^ FE + (M - 1,(51) 


9£ Sr] U e d£ 


Sr] j 


, 9 L 8E\ r (2^) 2n DE 

+ R 6 2 ^2 


(28) 


and for continuity: 


(2£) 2il || + |^ + (2?) 2R F^| + || log e 2^j = 0 (29) 

9uo 

In order to derive equation (29), it is necessary to obtain from equation (26) and the 

expression (from eqs. (22) and (24)) 


_ (2i|) n 9 

1 

CD 

& 

1 

\9XjJ v 9? 

(2^_ 


is also required. 

The continuity equation is written as shown to facilitate the computational procedure 
of Blottner (ref. 40) which is to be used herein. For h = — , equation (29) reduces to the 
form given by Blottner. 


Similar Solutions 

It is well-known from experimental investigations of turbulent boundary layers 
with dp e /dx = 0 that portions of the velocity profiles are at least approximately similar 
in terms of law -of -the -wall and velocity defect variables (ref. 41). As the Reynolds 
number approaches infinity, these profile shapes should approach exact similarity across 
the entire boundary layer. It is therefore appropriate to inquire as to the possibility of 
similar solutions to equations (27) to (29), since these equations should yield approxi- 
mately similar solutions if the Glushko formulations for the correlation terms are to be 
generally valid. 

To determine whether similar solutions exist, the question is posed; what distribu- 
tions of U e and n are required for F, V, E, and <p to be functions only of rj? 
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Inspection of equations (27) to (29) indicates that even with U e constant, no such exact 
solutions are possible because of the Rg factor in the M and D coefficients and in 
the last term (dissipation) of equation (28). This term may be written as 


2 n 


Dissipation 


_ c (2g) 1 DE _ 


R, 


r 


DE 

2.2 

V 6 <P 


since from equation (22), t ) g = \5 e &jv{2^) n . Substitution for D from equation (19) with 
R » 1 then gives the approximation 


E 3 / 2 U e 6 

Dissipation = C oik — — 

2 , v 
V 6 <P 

Thus, even if tj g is constant, the dissipation term does not satisfy similarity require- 
ments. However, if Rg is large and approximately constant over some small interval 
of interest and U e is constant, locally similar solutions would be possible. In general, 
exact similar solutions are not possible with the present variables £, 77 (eqs. ( 21 ) and 
(22)) and the transport coefficients M and D modeled according to equation (19). The 
existence of similar solutions is investigated further in appendix A where it is shown that 
if 6 increases linearly with Xj, similar solutions in terms of the variable x 2 /6 can 
be obtained for R » 1 and either U e constant or 

6 dU e „ 

— — - — = Constant 
U e dx t 

For Cj approximately constant, this condition is equivalent to that of Clauser (ref. 42) 
for equilibrium flows. It is also shown in appendix A that the ratio of a dissipation 
microscale to an integral scale should increase as some small power of the Reynolds 
number in order to obtain similar solutions of the flat-plate boundary layer. Actual solu- 
tions to the equations with similarity requirements included were not obtained in this 
investigation because numerical procedures would still be required, and little additional 
information could be expected beyond that just discussed. 

RESULTS AND DISCUSSION 
Computation Procedure 

The system of equations (27) to (29), along with auxiliary functions for M, D, 
and as defined in the previous section, is solved by a linear implicit finite- 

difference procedure. Complete details of the procedure including the finite-difference 
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expressions for the partial derivatives are given in appendix B. To increase the effi- 
ciency and accuracy of the computations, a variable grid size in the 7?- direction has been 
incorporated in the finite -difference expressions and used in some of the solutions. 
Details of the computer program including an operational flow diagram and program 
listing are given in appendix C. 


General Comments 

The effects of some modifications to the method of Glushko (ref. 17) on both mean 
and fluctuating flow properties are presented for flat-plate flow (dp e /dx = 0 ) and for one 
of Goldberg's experimental flows (ref. 43) with a large adverse pressure gradient. The 
principal modifications considered are to the l/b function (eq. (8)) and to the dissipation 
and diffusion terms (eqs. (12) and (13)). Limited results that indicate possible advan- 
tages to be gained by the use of similarity-type variables (eqs. (21) and (22)) with ii 
adjusted to give approximately constant boundary -layer thickness in the transformed 
coordinates are also mentioned. 

The data from the investigation by Goldberg (ref. 43) chosen as a test case should 
be a particularly severe test of any method (as indicated by Bradshaw and Ferriss, 
ref. 44) because the boundary layer was first driven nearly to separation by a large 
adverse pressure gradient and then allowed to relax toward a flat-plate flow by imposing 
a constant pressure. Another advantage of using reference 43 as a test case is that hot- 
wire measurements of turbulent shear and longitudinal turbulence intensity were made. 
Comparisons between these data and the theoretical results assist in the evaluation of the 
models used for the turbulent correlation terms. Also, the turbulence intensity data can 
be used to provide initial E profiles which are required in order to start the 
calculation. 

Bradshaw and Ferriss (ref. 44) and Nash (ref. 19) have applied their methods using 
the turbulent energy equation to this same adverse-pressure-gradient case of Goldberg 
and obtained reasonable agreement with experimental data for 9, H*, and Cj. Nash 
also obtained reasonable agreement between the calculated profiles of mean velocity and 
shear stress and the experimental results of Goldberg. Bradshaw and Ferriss (ref. 44) 
show that the computed values of Cj in the region of its minimum were very sensitive 
to the initial value of 6 used to start the calculation for this particular case. 

In all calculations by the present method, the skin friction has been computed from 
the correct limiting form evaluated at the wall as 
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hence from equation (22), 




= (2£)~ n 



and in finite -difference form 




E m,n=2 

A V=i 


(30a) 


(30b) 


In order to obtain valid results, the value of Ai7 n _i in equation (30b) must be less than 
the thickness of the region where iij varies approximately linearly with X£. 


Flat-Plate Flow 

The calculation was started at = 10 4 with input values of F and V from 

exact numerical solutions to the laminar Blasius flow. The input profile for the turbulent 
kinetic energy was taken (from ref. 17) as: 


E(S 0 *) " E o* 



(31a) 


where E Q * and 77* /rj e are specified constants, 
shown herein were computed with E Q * = 2.5 x 10 -4 
A77 = 0.05. Additional required inputs were 


Unless otherwise noted, the results 
* 

, = 0.4, ii = 0.5, and constant 

Ve 


U e = 100 ft/sec (30.48 m/s) 

77 e ,o = 4 - 9 5 

v = 1.58 x 10" 4 ft 2 /sec (l4.68 x 10" 6 m 2 /sec) 
e e = 1 x 10" 4 
e e ' = 3.4 x 10" 4 
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The Ai; step size was increased with increasing £ according to the relation: 

A| = f x 10 k (l.O x 10 k ^ k = 10 x 10 k ) (31b) 

so that the number of A£ steps required to traverse one cycle (on a log scale) of ij 
was the same for a given problem. Examples are shown for f = 0.1 and f = 0.5. 

The 1/ 6 function was obtained by linear interpolation from values given in table I. 
Three functions are given in table I: the one denoted as <pQ gg is based on the corre- 
sponding function as used by Glushko (ref. 17), and this <pQ gg function is used for the 
present calculation of flat-plate flow. The subscript 0.33 denotes the maximum value 
of l/6. 

The values of the constants in the transport functions M and D are: a = 0.2, 
k = 0.4, C = 3.93, R 0 = 110. As noted previously, these values and the <p function were 
adjusted by Glushko to give agreement with flat-plate flow. In the present report the same 
values, except where noted, are also applied to the adverse-pressure-gradient flow. 

Form factor, skin friction, and mean-velocity profiles.- The variation of the form 
factor H* with R x ^ is shown in figure 1. The value of H* is at first approximately 
constant at the initial value of 2.592, corresponding to the Blasius solution for laminar 
flow, and then H* abruptly decreases at some value of R x ^ depending on the step size 
(or f-factor), the magnitude of the input disturbance E 0 *, and modifications to the 
Glushko models of the correlation terms. This value of the Reynolds number where the 
flow proxies began to change u designated R^. Became of the behavior of H* as 
well as the mean-velocity profiles and other flow properties to be presented in subsequent 
sections, R^ ^ can be considered analogous to a transition Reynolds number.^ The 
two modifications used here were to the H(R) function (eq. (9b)) and to the diffusion term 
(eq. (13)). The reasons for investigating these modifications are discussed after their 
effects on H*, C^, and the mean velocity profiles are given. 

The modification to the H(R) function was to disregard equation (9b) in the outer 
part of the boundary layer and set H(R) = 1 when R/R q < 1.25. This modification had 
no effect (within the plotting accuracy) on the values of H*. 


^It must be emphasized, however, that no claim is made regarding the possibility of 
computing a realistic value of transition Reynolds number from the present method in its 
current stage of development. The term transition Reynolds number is used herein 
merely as a matter of convenience. In this connection, however, Donaldson has shown in 
a recent paper (ref. 45) that by retaining the complete equations for each component of 

the correlations u^'u. ' , some of the physical characteristics of transition can be com- 
puted. Donaldson has pointed out that his preliminary models for the correlation terms 
are crude and that improvement of the models should ultimately result in quantitative 
predictions of transition. 
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When the step size factor was increased to 0.5 (with H(R) = 1 for X 2/6 > 0.5), 
the value of R^ where H* first began to drop was increased slightly. It is apparent 
that this slight increase in R^ was a numerical effect caused by the larger step. 

Consequently, the following discussion of flat-plate results is limited to those solutions 
with f = 0 . 1 . 

Increasing the diffusion term by a factor of 3 increased the value of R x ^ by a 

factor of about 1—. Also when the input disturbance level E n * was reduced to 1 x 10“® 

3 4 

from 2.5 x 10"^, Rjj was further increased to about 8 x 10 . However, for all modi- 

fications shown here, the H* curves appeared to approach a value of approximately 1.4 
which is in agreement with the data of Wieghardt and Tillmann (ref. 46) considered typical 
of flat-plate flows. 

The effects of the modification to the diffusion term and the change in E 0 on Cf 
are shown in figure 2 and are of the same nature as the effects on H*. (The H(R) 
modification had a negligible effect on both Cf and H*.) That is, the "transition" 
Reynolds number is increased by the same factors when the larger diffusion term or the 
smaller E Q * are used but the final asymptotic values of Cf are in reasonable agree- 
ment with those for fully turbulent flow as obtained from Schlichting (ref. 39, p. 540) and 
the data of Wieghardt and Tillmann (ref. 46). According to the variation in Cf, when the 
boundary layer was fully turbulent at R^ ~4 x 10®, the value of H* was still 
decreasing at this Reynolds number as were the data of Wieghardt and Tillmann. (See 
fig. 1 .) 

The computed values of the mean velocity are plotted in conventional profile form in 
figure 3(a) and in law -of -the -wall and velocity defect coordinates (ref. 41) in figures 3(b) 
and 3(c), respectively. From figure 3(a) it is seen that the profiles develop from the 
laminar input profile at R x ^ = 1 x 10 ^ to turbulent type profiles for R x ^ 1 1 X 10 6 . At 
this Reynolds number the profile shapes have apparently settled out to the shape charac- 
teristic of turbulent boundary layers as indicated by the data of Wieghardt and Tillmann 
(ref. 46). In the transition region, that is, for 3 x 10^ < R x ^ < 2 x 10®, the profile shapes 
appear to be in qualitative agreement with mean profiles observed in transition flow as 
shown, for example, in reference 47. 

The effect of the H(R) modification (h(R) = 1 for X 2/6 > 0.5) is shown for the 
profile at R x ^ = 1 x 10®. This change in H(R) increases the mean velocity somewhat 
in the outer region of the boundary layer, with the result that the profile is fuller than 
when the original Glushko function for H(R) (eq. (9b)) is used. 

When the turbulent diffusion term is increased by a factor of 3, the agreement with 
data is considerably improved, again in the outer part of the boundary layer. These same 
effects are also shown in figure 3(c) (velocity defect region). The general agreement 
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with the law-of-the-wall correlations (fig. 3(b)) shows that the results in the near-wall 


region 


u*X2 


< 30 and the law-of-the-wall region are in good agreement with data. 


Discussion o f modific ations to H(R ) function and diffusion term.- The modification 
to the H(R) function consists of setting H(R) = 1 when R/R 0 < 1.25 in the outer part 
of the boundary layer rather than using equation (9b) in this region. (The H(kR) func- 
tion was modified in the same way, that is, H(«R) = 1 when kR/R 0 <1.25 in the outer 
part of the boundary layer.) The main effect is then to increase the eddy viscosity 
(eq. (9a)) for R < 137.5 (=1.25R 0 ); thus, the resulting fuller velocity profile shape as 
shown in figure 4 is to be expected. However, examination of the distribution of R 
across the boundary layer at R^ = 1 x 10® indicates that R > 137.5 for 
0.06 < X2/6 < 0.9. Hence, there was no increase in e in this range of X2/6 and yet 
the mean velocity was affected in the entire midregion of the boundary layer. This result 
merely shows that a local change in e can affect the computed mean velocities in the 
entire boundary layer, as it should for a nonsimilar numerical solution. 


The H(R) modification was thought to be more realistic than the original formu- 
lation because in the outer part of the boundary layer, any decrease in e should be taken 
care of by the boundary condition E e — 0 which, in turn, is caused by the intermittent 
properties of the flow rather than by a change in the basic relations for the turbulent cor- 
relations. In any event, since this modification to H(R) appeared to improve the veloc- 
ity profile shapes and had little effect on skin friction and form factor, the modification 
has been used in most of the remaining solutions presented herein. 


In regard to the modification of the diffusion term (3 times eq. (13)), a plot of the 

ratio of the turbulent shear stress to the turbulent velocity correlation, Tr^j2pe is shown 

in figure 4. Results computed from the theory illustrating the effects of the H(R) and 

diffusion-term modifications are compared with experimental data in this figure. The 

significance of this ratio and its tendency to be constant for widely different types of fully 

turbulent shear flows is discussed, for example, in references 9, 34, 48, and 49. The 

results obtained with terms in the turbulent kinetic energy equation identical to those of 

Glushko (ref. 17) were therefore considered to be unsatisfactory because of the large 

2 

increase in the shear ratio in the outer part of the boundary layer as shown in figure 4. 
The H(R) modification increased r,p in the outer part of the boundary layer so that the 


tion term 


^It was originally thought that the H(R) modification would increase the produc- 
9u-i\ 

in equation (2) and therefore cause an increase in e in the 


w 


0Xt 


outer part of the boundary layer. This increase in e should cause a reduction in 
Trp/2pe. Since the modification to H(R) was also applied in the same way to H(kR), the 
dissipation term was also increased. The net effect was only a slight increase in e for 
x 2 

■4 > 0 . 8 . 

0 
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resulting values of T,j,/2pe (fig. 4) were much too large for > 0.8. Since dF/8r] 

and E both approach zero near the outer edge of the boundary layer, the production and 
dissipation terms (the second and fourth terms on the right-hand side in eq. (28)) both 
become small there, and it is therefore apparent that another possible way to correct the 
T T /2pe ratio is to modify the diffusion term (the third term on the right-hand side of 
eq. (28)). ^ The simple expedient of increasing this term arbitrarily by a factor of three 
produced the desired result for the ratio Trj,j2pe as shown in figure 4. 

x 2/ 6 F=0.995 for 

several values of R x ^ and with the same modification to the diffusion term as noted 
previously. Results from two experimental investigations at R^ « 5 X 10 6 are shown 
for comparison. (The discrepancies between the data of refs. 32 and 50 are probably 
within the accuracy of the data.) When the terms in the turbulent kinetic energy equation 
are the same as those of reference 17, the changes in \[E profiles with Reynolds number 
duplicate those of reference 17 with some dependence on the peak turbulent intensity E Q * 
at the initial station = 4 x 10^. That is, when E 0 * is decreased, the peak values 
of \fE at subsequent stations in the transition region are also decreased (shown in 
fig. 5(b) with Glushko diffusion increased three times) but the profiles for Rj^ > 5 x 10^ 
are of the same shape and magnitude regardless of the value of E Q *. This computed 
behavior of the profiles is qualitatively in agreement with experimental observations 

in the transition region. (See ref. 47, for example.) Multiplying the turbulent diffusion 
term by three decreased \[e in most of the boundary layer except near the outer edge 
where \[E was at most doubled. These latter effects would be expected because of the 
gradient-type model (eq. (13)) used to formulate the diffusion term. 

Since the H(R) modification together with the arbitrary increase of three times 
the diffusion term improved the F and E profiles as well as the Tr^j2pe distribu- 
tion, the remaining discussion for the flat-plate case, as well as the adverse-pressure- 
gradient case, is limited to those results with both modifications included. It is empha- 
sized that these simple modifications appeared to improve the agreement between the 
calculations and experimental data of both the fluctuating and mean characteristics of the 
fully turbulent boundary layer. 

Bou ndary-lay er-thickness parameters .- The momentum thickness and the boundary- 
layer-thickness Reynolds numbers obtained from the solutions are plotted against R x _ ) 
in figure 6. The momentum thickness (fig. 6(a)) follows the laminar trend of |^R^~ 

3 

Other possibilities are to restrict the H(R) modification to the turbulent shear 
term only or to replace C in the dissipation term by a suitable function of X2/6. The 
diffusion term is also a candidate for modification because it is the only other term in 
equation (28) subject to semiempirical modeling. 


In figure 5 the computed values of {E are shown plotted against 
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until somewhat downstream of the values noted for the H* and Cj results. 

l,t 

Thereafter, 6 increases and appears to approach asymptotically the nominal variation 
of (ref. 39, p. 537) for turbulent boundary layers. The agreement with the 

data of Wieghardt and Tillmann (ref. 46) appears to be reasonable. 

The variation in Rg (fig. 6(b)), where 6 was taken as the value of X2 at which 
F = 0.995, shows the same trends as the 6 variation except that ds/dxj appears to 
increase more rapidly than dp/dxj and reaches a well-defined peak value in the neigh- 
borhood of R x ^ = 3 X 10^. This same type of variation in boundary -layer thickness within 

a transition region has been observed in experimental investigations such as that of ref- 
erence 47. These data show that the slope dS^dx^ is larger in the transition region 
than in the turbulent region to which it would eventually asymptote, that is, the slope 
d log e log e xj must approach the nominal value of 4/5 downstream of transition. 

This value of 4/5 is characteristic of fully turbulent boundary layers and appears to be an 
asymptotic limit for the present calculation. 

Balance of turb ulent kinetic energy .- The convection, production, diffusion and dis- 
sipation of turbulent kinetic energy corresponding to the first three terms, the fourth 
term, the fifth term, and the sixth term, respectively, in equation (28) are plotted against 
X 2/^F=0 995 * n fi S ure The scales used for the ordinate axes in this figure correspond 
to the dimensionless parameters used by Klebanoff (ref. 32). Computed values from the 
theory are shown in figures 7(a) and 7(b). The convection and diffusion terms are negligible 
in the region 0.1 < x 2 /& < 0.6. In the outer part of the boundary layer (X2/6 > 0.85^, these 
two terms are dominant over the other terms. In the inner part of the boundary layer for 
X2/5 < 0.01 (fig. 7(b)), the diffusion is of the same order as the production and dissipation 

terms and, as would be expected, the diffusion changes sign (proceeding outward from the 
wall) from a gain of energy very close to the wall to a loss of energy for 
0.001 < X2/6 < 0.025, back to a gain for 0.025 < ^ 2 /^ < 0-1- The production and dissipa- 
tion are essentially equal and opposite in sign for the region 0.1 < X2 fb < 0.6, and are 
still the largest terms in the near-wall region. Comparison of these results with those of 
Klebanoff (ref. 32) replotted in figures 7(c) and 7(d) indicates, in general, qualitative agree- 
ment except for the magnitude and trends of the diffusion for X2/5 < 0.5. The computed 
diffusion appears to be too small for X2/5 < 0.5 and also does not duplicate all trends of 
the experimental results. It should be noted, however, that there may be some error in 
Klebanoff's diffusion (which was not measured directly) in this region, since the integral 
of the diffusion from his data would not be zero, as it should be. Comparison of the com- 
puted diffusion (fig. 7(a)) with the measurements of Townsend, as quoted by Hinze (ref. 1, 
p. 498) shows that the trends in the computed values are in agreement with Townsend's 
data for Q.l < < 1.0. The agreement between the computed diffusion and the data 
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of Townsend could probably be further improved by multiplying the original diffusion 
term by a simple function of X£ rather than by a constant. 

Nonequilibrium Boundary Layer 

Several investigations (refs. 36, 43, and 51, for example) have been made recently 
to determine the response of a turbulent boundary layer to changes in the external pres- 
sure gradient. In the investigation of Goldberg (ref. 43), detailed hot-wire measurements 
of turbulent shear stress and longitudinal turbulence intensities, as well as measurements 
of velocity profiles and wall shear stress, were made for six different pressure distribu- 
tions. The test boundary layers of reference 43 were established on a 10-inch-diameter 
cylinder which was alined parallel to the free stream. The flow was carefully checked for 
axisymmetry from lateral variations of wall static pressures, wall shear stress measured 
with sublayer fences, and the separation line as indicated by tufts and the sublayer fences. 

Since the maximum boundary-layer thickness was about 1.5 inches (0.04 m), the 
effects of the axisymmetric geometry on the boundary -layer characteristics must be con- 
sidered. Goldberg found that the maximum difference between values of 6* and 6 from 
axisymmetric and two-dimensional formulas was about 10 percent whereas the maximum 
difference for H* was only about 2 percent. Differences in skin friction and fluctuating 
properties would probably not be any larger. It is considered likely that differences of 
this type are smaller than three-dimensional effects in a conventional two-dimensional 
test channel, particularly if the boundary layer is near or at separation. The computer 
program described herein is at present limited to two-dimensional flows. The additional 
complexity of accounting for the axisymmetric geometry was not considered to be worth- 
while in view of the preliminary nature of the models used for the various turbulent 
correlations. 

It was concluded that from among the available investigations, the data by Goldberg 
would be the most useful as test cases and his pressure distribution number 3 (ref. 43) 
was selected as the most severe test case since, as mentioned previously, the boundary 
layer was driven almost to separation by a severe adverse pressure gradient over the 
first 16 inches (0.40 m) of the run and then allowed to relax toward a flat-plate boundary 
layer by imposing a constant pressure run for the next 24 inches (0.61 m). 

The distribution of external velocity and its derivative with respect to Xj as used 
to obtain the present results are shown in figure 8. There was some uncertainty in 
reading the small graphs published in reference 43. Also, this case appears to be unusu- 
ally sensitive to initial conditions as shown by Bradshaw and Ferriss (ref. 44). Conse- 
quently, all input quantities as used in the present calculations are listed in tables n 
and III. Second-order interpolation was used between these tabulated values to supply 
sufficient detail for the calculation. Table II lists the values of F and E used at the 
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input station corresponding to = 4 inches (0.10 m) (presumably the distance from the 
nose of the 10-inch-diameter (0.25 m) test cylinder). Table III includes U e , dU e /dxj, 
dUe/d^, and £ as functions of xj_. Two different tabulations, designated as velocity 
distributions a and b, are given in table in. These alternate velocity distributions are 
plotted in figure 8(a) and the differences between them are believed to be within the reading 
accuracy of the original figures in reference 43 and probably within the experimental 
errors of the original data. The resulting derivatives of the two velocity curves are con- 
siderably different as shown in figure 8(b) where dU e /dxj used in the calculations is 
plotted against xj. These differences can become important when the boundary layer 
approaches separation. 

The initial | 0 (8 x 10^) for the input station (x^ = 4 inches (0.10 m)) was computed 
by assuming flat-plate flow at U e = 85 ft/sec (25.9 m/s) over the upstream test cylinder 
length taken as 1.5 feet (0.46 m). The solutions should not be sensitive to the values of 
£ 0 if all other input conditions are held fixed. This insensitivity to was verified by 
obtaining a repeat solution with | 0 = 3 x 106, and all results were essentially identical to 
those obtained with ^ 0 = 8x 10^. The initial velocity profile was taken directly from the 
data plot of reference 43 at xj = 4 inches (0.10 m) and the initial E profile was taken 
from the measured longitudinal intensity at the same station as 


E = 



+ Ur 


»2 


2U f 




(32) 


This expression is based on the data of reference 32 where all these components of the 
fluctuating velocity were measured. These data are plotted in figure 5(a) which shows 
that the approximate expression (32) is accurate for 0.1 < X2/5 < 0.8. A constant value 
of A77, corresponding to K = 1, was used for most of the solutions reported herein. 

Repeat solutions with variable A 77, corresponding to K = 1.02 (see appendix B) gave 
essentially identical results and reduced the number of A77 steps to about 1/4 of the 
number required for K = 1.00. 

Relation between skin friction and \fE profiles.- The calculated variation in skin 

friction is compared with the experimental data in figure 9. The shaded band represents 
the spread in the experimental skin-friction data obtained by four methods (ref. 43): 

(1) sublayer fence, (2) Preston tube, and from measured velocity profiles using (3) Clauser’s 
method (ref. 42) and (4) the Ludwieg-Tillmann relation (ref. 52). Theoretical results are 
presented for both velocity distributions of figure 8, three values of C (dissipation con- 
stant), and three 1/ 6 functions. The (p q gg function is the same as that used for the 

flat-plate calculations given in table I and the other two functions as used herein are also 
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given in table I. Linear interpolation between the tabulated values was used in the solu- 
tions. Prandtl's mixing length relation is thereby recovered in the law-of-the-wall region 
(for example, X 2/6 < 0 . 2 ) where production and turbulent dissipation are approximately 
equal. That is, by equating production and dissipation 


9uj 9uf* 9uj' 
T T 9x 2 ^ 0x j 9x j 


(33) 


and by eliminating e between equations (7) , (9a) , and (12) with R » 1 (H(R) and 
H(kR) = 1), there is obtained 


II « « > 9 . /^ u ! \ 
P \[kC \ 9 x 2 / 

Then with l = y, a = 0.2, k = 0.4, C = 3.93 


(34a) 


T ~f - (0.4y ) 2 



(34b) 


which corresponds to the Prandtl mixing-length relation for turbulent boundary layers. 

In the following discussion, the effect of velocity distribution and C are considered first. 

For velocity distribution a, C = 3.93, and 1 / 6 = <£g 33 , the agreement with data is 
good (fig. 9) except in the region of the minimum Cf where the present method over- 
predicts Cf by as much as 100 percent. The agreement is somewhat better when the 
alternate velocity distribution b was used. These different values for Cf in this region 
indicate the sensitivity of the results to the imposed velocity distribution. 

In order to determine the relative importance of the dissipation term (eq. (12)) for 
this particular type of flow, additional solutions with C = 5 and C = 6 for pressure 
distribution b and <p q 33 were obtained. An increase in C increases the dissipation 

and reduces the skin friction by an almost constant amount over the entire test region and 
gives improved agreement with data near the minimum Cf region, at the expense of 
poorer agreement elsewhere. The turbulent kinetic energy equation can be paraphrased 
as 

^ = Production + Diffusion - Dissipation (35) 


where ^ is the derivative following the particle. It is evident that an increase in 
dissipation should decrease e, and this decrease, in turn, should decrease the turbulent 
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shear which from equation (10) is approximately 


t T 



( 36 ) 


for R > R 0 . Since the skin friction depends directly on the magnitude of in the wall 
region, the decrease in Cf due to the increase in C appears to be reasonable, although 
it is somewhat surprising to find the almost linear relation between the magnitudes of the 
dissipation and Cf. 


In order to arrive at an explanation for the improved agreement between the com- 
puted and experimental skin friction in the minimum Cf region, as caused by larger 
values of C, profiles of \Je are shown in figure 10. The theoretical values in this fig- 
ure were computed with the diffusion term taken as three times the Glushko diffusion 
(3 x eq. (13)), velocity distributions a and b, different values of C, and different l / 6 
functions. The first thing to note is that all these various modifications had only minor 
effects on the magnitude and distribution of at a given xj station. An increase in 


C 


does reduce 



as it should according to equation (35), but the best agreement 


with the data is generally obtained with C = 3.93. (The effect of the change in the l/& 
function is considered in a subsequent section.) The input profile of \/e (fig. 10(a)) is 
similar to the flat-plate profiles of figure 5. Comparison of the \Je profiles at 
increasing Xj stations as shown in figures 10(a) and 10(b) shows that as the minimum 
Cf region (near xj = 20 inches) is approached, the peak in both the computed and mea- 
sured ]/E profiles moves away from the wall and increases in magnitude. Consequently, 
there is a corresponding increase in the average turbulent energy (e) av across the entire 
boundary layer as x^ increases. That this increase in (e) av can be associated with a 
decrease in a dissipation length scale can be seen by noting that for large R, the 
dissipation (eq. (12)) can be written approximately as 


0X i 


8 ^ 

8 Xj 


vC 


1 +£(kR) 

- ^ J Zd 


~ Ko^e) 3 / 2 l -Rj 
2 Id 2 


(37) 


where ^ is defined here as a microscale of the turbulence (ref. 1, p. 37). Then if the 
quantity C jl^ 1 is considered as an adjustable function but with C fixed as a universal 

constant, the larger values of C as used for the solutions of figure 9 should be considered 
as equivalent to corresponding decreases in the square of the dissipation scale The 

quantity l in equation (37) is considered to be the same mean scale (analogous to an 
integral scale) as defined by equation (8) and as used by Glushko (ref. 17). 
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The relation between l and Z d for isotropic turbulence (ref. 1, p. 185) may be 
written as 



(38) 


where C T is a constant. This relation is seen to be in agreement with the results of 
figure 9 as related to the change in \/E profiles with xj (fig. 10) since as (e) av 
increased, l d should decrease which, according to the preceding reasoning, accounts for 
an increase in dissipation and the reduction in Cf, as computed. The reduction in Cf 
follows from the approximate equivalence of production and dissipation (eq. (33)) but with 
the dissipation given by relation (37). Relation (34a) then becomes 


Jt ~a_ J du l\ id 
p \Jk \Bx 2 ) \fc 

and substitution of the expression (38) into this relation gives 


II ^ 

P )/« Ue 




which indicates that a functional relation for Z d of the type given by equation (38) would 
improve the agreement between theory and data over the entire test length for this case 

because of the way (|/e) varies with xj. This explanation depends on the assumption 
that the mean scale l is not affected by (e) av ; that is, 1 / 6 is assumed to be a fixed 
function of x 2 /5- 

Boundary- layer thickness par ameters and mean velocity profiles .- The computed 
values of 6 and 6 are compared in figure 11 with the experimental data. In figure 11 
and all remaining figures, the theoretical values of 6 are taken at the point where 
F = 0.995, and the experimental values of 6 were taken directly from figure 20 of ref- 
erence 43. The differences between the theoretical and experimental values of 6 are 
small. However, these small differences can lead to large effects on Cf since the 
momentum integral equation 


dd _6_ dUe 

dxf Ue dxj 


(H + 2) = 


(39) 
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shows that when Cf approaches zero in a large adverse pressure gradient (negative 
dUe/dxi), the skin friction equals the difference between two large terms, and any small 
changes in either of these terms can lead to large changes in Cf. The difference 
between the theoretical and experimental values of 6 are also small, as shown in fig- 
ure 11(b). The use of velocity distribution a, C = 3.93, and 1/ 6 = 0q 25 generally 
gives the best agreement with the data except for the 6 distributions when xj = 20 inches 
(0.51m) where the 0 q. 33 solution (with C = 3.93 and velocity distribution a) gives the 
best agreement. However, it is obvious that a comparison of 6 and 6 alone are not 
sufficient to judge the accuracy of a method when Cf approaches 0. 

The shape parameter H* is a more sensitive indicator of the accuracy of a method 
as shown by figure 12 where the effects of the two velocity distributions, the values of C, 
and the l/5 functions are shown. The velocity distribution b and the largest value of 
C give the best agreement with the experimental data for 1/5 = 00.33- However, the use 
of the other two l / 5 functions gives the best overall agreement with the data, and the 
results bracket the data in the vicinity of the peak H*. 

The reason for the better agreement of H* with data for l/5 = 0 q 2 q and 0 q 25 
is apparent from figure 13 where the computed velocity profiles at three xj stations are 
compared with the data. (Here again the subscripts denote the maximum values of the l/ 5 
functions; these functions are given in table I.) The agreement between theoretical and 
experimental velocity profiles for all 0 functions is fair at xj = 12 inches (0.30 m) 

(fig. 13(a)) where dU e /dXf, according to figure 8 (b), is near its peak value. However, at 
xj = 20 inches (0.51 m) (fig. 13(b)) where dUe/dXf is approaching zero, the agreement 
is poor for l/5 = 0q 33 for both velocity distributions a and b and for all values of 
C. When 0q 20 and 00 25 are use d, the theoretical results bracket the data for this 
station of x^ = 20 inches (0.51 m), which is apparently a critical region where the cal- 
culation is very sensitive to the 1/ 5 function. Near the end of the constant pressure run 
at Xj = 36 inches (0.91 m) (fig. 13(c)), the good agreement between the theory (with 
l/5 = 0 q 33 ) and data may be fortuitous. 

Discussion of tur bulence scale functions.- It has already been noted from figure 10 
that \JE is relatively insensitive to changes in free-stream velocity distribution, C, and 
the 1/5 function. Figure 13 shows that the mean velocity profiles at a given Xj station 
are insensitive to the two velocity distributions and the values of C considered, but can 
be affected considerably by a change in l/5. Since the turbulent shear, production, and 
dissipation are the dominant terms in the equations, it is apparent from the preceding dis- 
cussion and from the form of these terms (see eqs. (36) and (37)) that the only way to 

^Note that in the present method Cf was not computed from equation (39), but from 
equation (30) which probably requires more stringent accuracy criteria in the calculation 
than equation (39) would. 
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change the mean velocity profiles, and hence the H* values to any appreciable extent is 
to modify the turbulence scale function 1/ 8. In the previous section, a change in C was 
related to a change in a dissipation microscale but from equation (34), a change in C 
is equivalent to a change in the mixing length constant of Prandtl. Since the mixing length 
relation of Prandtl is known to apply even in an adverse pressure gradient in the law-of- 
the-wall region, the 1/ 8 function was not changed for X 2 /S < 0.2. Goldberg's experi- 
mental values of mixing length (ref. 43) indicate that l/ 8 should be decreased in the 
region of X 2/6 > 0.2. Since this decrease would reduce the turbulent shear and thereby 
result in a more linear velocity profile, the 1/ 6 functions were changed in the manner 
shown in figure 14(a). (See also table I.) 

The corresponding values of the eddy viscosity function e/p at Xj = 20 inches 
(0.51 m) are shown in figure 14(b). As would be expected, the e/p function does not 
depend on the two free-stream velocity distributions or the values of C used, since ]fE 
was independent of these parameters. However, the distribution and magnitude of e/p 
are directly dependent on 1/8 and the theoretical curve for l/8 = 4>q 20 is closer to the 
experimental values of reference 43. This direct dependence of e/p on l/8 is, of 
course, the reason for the marked effect of the 1/8 function on both Cf (fig. 9) and 
uj/Ug (fig. 13). 

The success of these modifications to l/8 indicates that further minor adjustments 
to these functions should give further improvement in agreement with experimental data. 

In particular, the trends in \/e (fig. 10) for this problem indicate the magnitude of the 
change in l/8 should probably depend on the level of (e) av . Spalding (ref. 53) has 
applied to free shear flows a differential equation for the mean scale of turbulence. This 
equation was based on Rotta's hypothesis (ref. 6 , part II) and relates the scale of turbu- 
lence to (e) av and r,p. It is possible that the use of a similar relation may improve the 
predictions of the present method. 

Comparisons with empirical co rrelations of mean velocity .- Comparisons between 
calculated and experimental values of mean velocity are also shown in figure 15 in law- 
of-the-wall coordinates. At x-f = 12 inches (0.30 m) (fig. 15(a)) both the data and theory 
deviate from the logarithmic law of the wall for large values of u* as is character- 

istic of boundary layers with adverse pressure gradients. In the minimum Cf region 
(xj = 20 inches (0.51 m), fig. 15(b)), the data and the theory for 0g 25 still follow the 

u* x o 

law of the wall up to — — ~ 70. At the downstream station of Xf = 36 inches (0.91 m) 
(fig. 15(c)), after a constant pressure run of about 12 inches (0.30 m), the logarithmic 

u*X 2 

variation is followed by both data and theory (with <£o.25) f° r 25 < — — < 400. It can 
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be concluded, therefore, that the mean velocity profiles have completely adjusted to the 
zero-pressure -gradient condition. 

Fluctuating flow quantities .- The total turbulent intensity profiles in the form of l/E 
are shown in figure 10 and have been discussed previously. Experimental values were 
obtained from the measured longitudinal intensities by the use of expression (32). The 
agreement in both trends and magnitude is generally good over the portion of the boundary 
layer where the approximate relation (32) is applicable. 

The variation in the ratio of turbulent shear to twice the turbulent kinetic energy 
across the boundary layer at three stations is shown in figure 16. As mentioned previ- 
ously, this ratio r j,j2pe tends to have an approximately universally constant value in 
the intermediate region of the boundary layer. By the use of equation (32), the ratio 
r-p/2pe was obtained from the data of Goldberg (ref. 43) and Bradshaw and Ferriss for 
a relaxing boundary layer (ref. 36). These values are plotted for comparison with the 
theoretical results in figure 16. The stations of x-^ = 47, 59, and 95 inches (1.09, 1.5, 
and 2.41 m) should be comparable to Goldberg's stations x^ = 14, 20, and 36 inches 
(0.35, 0.51, and 0.91 m) (figs. 16(a), 16(b), and 16(c)), respectively, in terms of the rela- 
tive distributions of external velocity and velocity gradient in the two investigations. The 
theoretical results are generally in reasonable agreement with the data. (Note that the 
best agreement with Goldberg's data was generally obtained with (1/ 5) = <pQ 2 q or <£o.25-) 
It is seen that the ratio Trp/2pe is roughly the same at all stations and for the two sets 
of data shown except at the critical station of xj = 20 inches (0.51 m) (fig. 16(b)). Note 
that at xj = 36 inches (0.91 m) (fig. 16(c)), where the mean velocity profile and H* 
have relaxed to flat-plate values (figs. 15(c) and 12), the experimental and theoretical val- 
ues of r i^pe are still somewhat below that typical of flat -plate flow. 

The distributions across the boundary layer of the four terms in the turbulent kinetic 
energy equation at four x^ stations are shown in figure 17 for velocity distribution a, 

C = 3.93, and 1/6 = 00.25- The scales used in the figure correspond to the dimensionless 
parameters used in equation (28). At xj = 8 and 12 inches (0.20 and 0.30 m) (figs. 17(a) 
and 17(b)), the distributions are similar to those of the equilibrium boundary layer of 
Bradshaw for a = -0.15 (fig. 14(c) of ref. 34). At xj = 17.4 inches (0.44 m) (fig. 17(c)) 
the trends are more nearly like those of Bradshaw for a = -0.255 (fig. 14(d) of ref. 34) in 
that the relative magnitude of the diffusion is increased in the central part of the boundary 
layer. The distribution of production and convection terms at x^ = 20 inches (0.51 m) 
and 36 inches (0.91 m) (figs. 17(c) and 17(d)) are qualitatively similar to the corresponding 
results for the relaxing boundary layer of Bradshaw and Ferriss at xj = 53 (1.35 m) and 
83 inches (2.11 m) (ref. 36). 

Comparison of the results shown in figure 17 with corresponding results for velocity 
distribution b and other values of C indicates that the trends and magnitudes of the 
various terms are not greatly affected by these changes in external velocity distribution 
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and C except near the wall. The location and magnitude of the peak values (nearest the 
wall) of the terms are shown in the following table: 


Velocity 

distribution 

1/5 

C 

x 2 /e 

Production 

x 2 /S 

Dissipation 

x 2 /& 

Diffusion 

x 2 /6 

Convection 





X 1 = 

8 inches (0.20 m) 





a 

^0.33 

3.93 

0.021 

1.37 

0.005 

-2.55 

0.005 

2.05 

0.005 

-0.445 

b 

^0.33 

3.93 

.021 

1.4 

.005 

-3.86 

.005 

2.85 

.005 

.919 

b 

^0.33 

5 

.026 

1.25 

.016 

-.833 

.0106 

1.0 

.0106 

-.597 

b 

^0.33 

6 

.026 

1.139 

.0107. 

-.785 

.0107 

.764 

.016 

-.224 





X 1 = 

20 inches (0.51 m) 





a 

^0.33 

3.93 

0.014 

0.533 

0.0028 

-1.08 

0.0028 

0.919 

0.0028 

0.154 

b 

^0.33 

3.93 

.014 

.483 

.008 

-.376 

.0054 

.534 

.0054 

-.352 

b 

^0.33 

5 

.017 

.384 

.0028 

-1,44 

.0028 

1.18 

.0055 

.258 

b 

^0.33 

6 

.017 

.315 

.0028 

-1.0 

.0028 

.76 

.0028 

.244 





X 1 = 

20 inches (0.51 ra) 





a 

^0.20 

3.93 

0.019 

0.139 

0.0028 

-0.778 

0.0028 

0.654 

0.0058 

0.137 

a 

^0.25 

3.93 

.017 

.293 

.0028 

-1.27 

.0028 

1.07 

.0058 

.246 


Comparison of these peak values shows that the dissipation is not always increased when 
C is increased. This result is presumably an indication of the nonlinear nature of the 
problem. It is also of interest to note that the location of the peaks in all terms is usually 
closer to the wall (in terms of X2/6) at x-^ = 20 inches (0.51 m) than at xj - 8 inches 
(0.20 m). Another item worth noting from the table is that when the velocity distribution 
and C are the same but 1 / 6 is different, the magnitude (except for the dissipation and 
diffusion) and location of the peaks are appreciably different, in spite of the fact that the 
l / 6 function was not changed in the wall region. (See fig. 14(a).) This result is, of 
course, an important confirmation of the experimental fact that the upstream history of 
the flow in this kind of a boundary layer has a large effect on local conditions. 

Other adverse-press ure-gradient cases .- Any method for predicting the development 
of turbulent boundary layers cannot be considered satisfactory until the results for several 
different types of flows with a variety of initial and boundary conditions are compared with 
experimental data. Some comparisons of this type for the modified Glushko method have 
been reported by the present authors in reference 54. The method was applied to several 
cases and reasonable agreement in H*, Cf, and Rg was generally obtained. 
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CONCLUDING REMARKS 


The equations for the incompressible turbulent boundary layer with constant fluid 
properties have been solved by a numerical procedure in similarity type coordinates. A 
parameter in the definition of these coordinates could be adjusted to keep the turbulent 
boundary -layer thickness in the transformed coordinate system approximately constant; 
this procedure gives a nearly constant number of computing steps across the boundary 
layer. The conservation equations for mass, mean momentum, and turbulent kinetic 
energy are solved simultaneously by a linearized, implicit, finite-difference procedure 
wherein all boundary conditions at the surface and outside edge of the boundary layer are 
satisfied directly. 

Mathematical models of the turbulent production, dissipation, and diffusion terms 
developed by Glushko (Bull. Acad. Sci. USSR, Mech. Ser., no. 4, 1965) for flat-plate turbu- 
lent flow have been modified and applied also to a nonequilibrium turbulent boundary layer 
subjected initially to a large adverse pressure gradient which is followed by a constant- 
pressure region. Comparisons of calculated values for both mean and fluctuating flow 
properties with experimental measurements in this nonequilibrium boundary layer as well 
as the flat-plate boundary layer have indicated generally good agreement. 

For the flat-plate calculation, the laminar Blasius velocity profile and arbitrary 
small disturbance-type profiles for the turbulent kinetic energy were used as initial con- 
ditions at a Reynolds number of 10^. As the calculation proceeded, little change in the 
mean profiles of velocity and turbulent kinetic energy were noted until at some downstream 
station, depending on the level of the input disturbance and a modification to the turbulence 
diffusion term, rather abrupt changes began. The subsequent mean velocity profiles, tur- 
bulent kinetic energy profiles, and skin friction were qualitatively similar to those 
observed experimentally in the transition region between laminar and fully turbulent flow. 
The dependence on the level of the input turbulent kinetic energy of the Reynolds number 
at which these changes in the mean profiles were first obtained in the calculation was 
shown by Glushko. It is shown herein that modifications to the models of the turbulent 
terms also affected this Reynolds number. The value of this Reynolds number was an 
order of magnitude lower than minimum experimental transition Reynolds numbers and 
cannot yet be related to the physical phenomenon of transition. Nevertheless, the indica- 
tions are that further development of the method may ultimately yield a useful technique 
for the numerical investigation of transition. 

It was found that when the turbulent diffusion term of Glushko was multiplied by an 
arbitrary factor of three, the agreement with experimental values of mean velocity and the 
ratio of turbulent shear to turbulent energy was improved in the outer portion of the fully 
turbulent boundary layer for both the flat-plate and nonequilibrium flows. It was also 
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found that when the dissipation term was changed, the skin friction was reduced by almost 
constant increments that depended directly on reductions in the square of the dissipation 
scale. Analysis of these results indicated that the microscale (dissipation scale) may be 
related to the integral scale of turbulence in about the same way as that for isotropic tur- 
bulence. The best overall agreement with measured values of skin friction, form factor, 
mean velocity profiles, and fluctuating properties was obtained by reducing the value of 
the turbulence scale from a peak of 0.33 to 0.20 or 0.25 of the boundary -layer thickness 
in the midregion of the boundary layer. The linear relation of turbulence scale with dis- 
tance from the wall, in accordance with Prandtl's mixing length theory and as used by 
Glushko, was retained in the law -of -the -wall region. It is concluded that simple modi- 
fications to the turbulence scale function and to the turbulent correlation terms as modeled 
by Glushko result in accurate calculations of mean and fluctuating characteristics of tur- 
bulent boundary layers with arbitrary boundary conditions. 

The existence of similar solutions to the equations with the turbulence correlation 
terms modeled according to Glushko was investigated briefly. It is shown that exact 
similar solutions are not possible in the coordinates used in the numerical calculation. 
However, approximate similar solutions exist for other coordinate systems. In one of 
these systems, the normal coordinate is simply the physical distance from the surface 
divided by the local boundary-layer thickness which is then required to increase linearly 
with distance along the surface in order to obtain similar solutions. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Langley Station, Hampton, Va., August 1, 1968, 

124-07-01-32-23. 
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APPENDIX A 


REQUIREMENTS FOR SIMILAR SOLUTIONS 


Analyses of results from many experimental investigations have shown that the 

mean velocity profiles in the incompressible turbulent boundary layer, with dp /dx-i = 0, 

u 1 x„u* u 1 -U x 2 ' 

are similar in terms of — ± — — and — — in the inner and outer parts of the 

u* v u * o 

boundary layer, respectively (ref. 41). For dp e /dx^ ^ 0, Clauser (ref. 42) has 
velocity profiles are self-preserving or approximately similar if the parameter 

is a constant. As the Reynolds number approaches infinity, all velocity and length scales 
except 6 and U e should become independent of xj; therefore, even with an imposed 
pressure gradient, the entire mean velocity profile would presumably depend only on x 2 / 6. 
For a limited range of conditions or within certain regions of the boundary layer, similar 
solutions to the turbulent boundary- layer equations are therefore possible as shown, for 
example, in references 9, 14, 20, 55, and 56. The purpose of this appendix is to discuss 
briefly the conditions required for the existence of similar solutions to the equations of 
mean momentum, turbulent kinetic energy, and mass continuity in the form of equa- 
tions (16) to (18). It is hoped that the discussion will assist in the evaluation of the models 
proposed by Glushko for the shear, dissipation, and diffusion of turbulence. 

For this purpose, it is convenient to introduce general coordinates defined as 


shown that 
/5*\dp e 
\ T wjdx^ 



y ( x 1’ x 2) = PQ x 2 


(Al) 


where P and Q are as yet unspecified functions of x-^. It is necessary at the outset 
to exclude the viscous sublayer region in order to obtain similarity. Hence, the restric- 
tion R » 1 is imposed and the coefficient functions M and D (eq. (19)) become 


M 


a 


v ' 6 


D « GIK 


Up 6 


IeI 


(A2) 


The transformation of equations (16), (17), and (18) to the variables X and Y then gives 
for momentum: 


F SF T7 - 9F _ 1 - F 2 dU e 
Q 3X + V 3Y Q Ue dX 


+ oPQ6 — 
8Y 



(A3) 
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for energy: 


F 9E 8E 
Q ax + V 9Y 


-2 


FE dU e 
QU e dX 



+ 



Caw 


E 3 ^ 2 

PQ? d 


and for continuity: 


l_aF 9V F_d_/ Ue\ 
Q ax ay Q dX^ s epQ/ 


= 0 


(A4) 


(A 5) 


where 


F ay u 2 
PQ Bxj + U e 


(A6) 


From equations (A3), (A4), and (A5), it can be seen that the profile functions F, E, and 
V may be functions only of Y if the following conditions are satisfied: 



1 

QU< 

dU € 

jdX 

3 1 

PQU e 

dU e 

dx-^ 

= c i 

(A7) 

1 d 

(log e 

u e \ 

1 d 

1 OP’ 


(A8) 

Q dX 

PQ/ 

PQ dX! 

l lu &e 

! pq; 2 




c 3 



(A9) 



l 

6 



(A 10) 


where the C values are constants. The transformed variable Y then becomes simply 

Y = C 3 ^ (All) 


and equation (A10) is the function cp as assumed by Glushko. The requirements on U e 
and 6 are obtained by combining equations (A7), (A8), and (A9) to give 


U e dxj 


= Ci 


(A12) 


^ = C 2 - Cl (A13) 
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Hence, 5 must vary linearly with in agreement with the results of reference 56 
(pp. 81-83). (See also ref. 7, pp. 105-109.) By means of Bernoulli’s equation, condi- 
tion (A12) may be written as 


5 dp _ ^ 

^ 7 ^' 1 


which is essentially the same as Clauser's condition for equilibrium boundary layers 
(ref. 42, pp. 49-50) 

— — — = Constant 
pU e 2 C f dx 

provided that 6*/5Cf is approximately constant as would be the case for equilibrium 
layers. Consequently, the requirement of equation (A12) would be approximately satis- 
fied by equilibrium boundary layers (including the zero-pressure-gradient case) but the 
requirement of a linear variation of 6 with x^ would not be satisfied in general. 

Another class of similar solutions without this restriction on 6 can be derived by 
dividing equations (A3) to (A5) by P. (Note that in the previous class of solutions where 
Y varies as x 2 /S, P is arbitrary.) The transformed normal velocity is then defined by 

V(Y)=-F-^ + -^_ 

P 2q axj PU e 


and the remaining requirements for similar solution are 


dU c 


P 2 QU e d*! 


= Ci 


(A 14) 


2n dl ( l0g e P q) C 2’ 


1 

P 2 Q 


(A 15) 



(A 16) 


where all C quantities are constants. The function 1/6 in equations (A3) and (A4) is 
again required to be a function only of Y 

(A 17) 

However, the dissipation scale l ^ in the last term of equation (A4) must now satisfy the 
relation 
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Zd = C4 

l p2 


where equations (A16) and (A17) have been used. 

To find the P function, equations (A14) to (A16) are combined to give 


8 u u _ n e 

P dx -P ^ 


d 6 


where C 5 = ^ 2 ' - C-j^Cg'. Integration of this expression yields 



6 dXj 


(A18) 


(A 19) 


and the ratio of the dissipation scale to the transport scale (eq. (A18)) becomes 


^d _ ^6 f x l 
1 6 2 ^0 


6 dx^ 


(A20) 


where Xj = 0 represents the effective origin of the turbulent boundary layer. The con- 
ditions given by equations (A16), (A19), and (A20) are more general than those of equa- 
tions (A12) and (A13) since no restriction on 6 is required. 

The P and 1$ functions given by equations (A19) and (A20) could presumably be 
determined by numerical iteration techniques for any boundary layer; however, equa- 
tions (A14), (A16), and (A19) give the additional requirement that 

C’6 dx. 

d(log„ UJ = — — 

6 dxj 


\ c 


f 

J 0 


as 


For 



and over a limited range of x^, the boundary -layer thickness varies 



(A21) 


Substitution of this expression into equation (A20) yields 



(A22) 


Thus, if n = 0.8, which is a realistic value, the dissipation length scale would increase 
as the 0.2 power of the Reynolds number in order to obtain similar solutions for a 
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flat-plate boundary layer with 5 increasing according to relation (A21). It is of interest 
to compare the ratio l^fl (ratio of dissipation microscale to the transport, or integral, 
scale) from equation (A22) to the corresponding ratio for isotropic turbulence as given by 
Hinze (ref. 1, p. 185, and eq. (57)). This result for isotropic turbulence may be written 
as 


l 


d 

l 


C C 



(A23) 


or by the use of relation (A21) 



(A24) 


Equation (A24) indicates that if l / 6 and E are similar profiles, and if the flat-plate 
dissipation varied in the same way as isotropic dissipation, l^jl would decrease with 
increasing Reynolds number. The isotropic expression (A23) appears to be in qualitative 
agreement with requirements for in nonequilibrium boundary layers (see discussion 
concerning eq. (38)); however, by comparison with equation (A22), it can be concluded that 
the isotropic expression would not be applicable to a flat-plate equilibrium boundary 
layer. 

Other forms for P and Q that would yield similar solutions can be derived by 
multiplying equations (A3) to (A5) by P m Q 1 where m and i are constants. Since 
these other forms appear to give results of limited interest or applicability, they are not 
discussed. 


42 


I 



I 


APPENDIX B 

NUMERICAL COMPUTATION PROCEDURE 

The system of equations (27) to (29), along with auxiliary functions for M, D, 
and n^x-jJ as defined in the section "Transformation to Similarity Coordinates" is 
solved by a linear implicit finite -difference procedure. This procedure combines cer- 
tain aspects of the methods given in references 17 and 40. 

Grid Notation 

The boundary -layer region is divided into finite grids of A£ width and A r\ height 
as shown in the sketch, where the grid or nodal -point notation is also illustrated. To 
increase the efficiency and accuracy of the finite -difference procedure, a variable grid 
size in the 77 -direction is included. The grid size in the ^-direction may also vary. 
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The variation in the A 77 grid size is controlled by the parameter K where, as 
illustrated, 


K = 


A T] 


n 


Ar ?n-1 


(Bl) 


When K is a constant, the successive Ar] n values form a geometric progression; 
hence, 


Ar, n = K 11 ' 1 A Vl (B2) 

The total number of steps across the boundary layer is N - 1; thus, the size of the last 
step is 


A77 


N-l 


= K N ‘ 2 Ar? 1 


(B3) 


The maximum thickness of the boundary layer r) e is given by 

K N-l 1 

% = — f- 4 ”i 


(B4) 


Thus if ? 7 e , Ar;^, and the number of steps N-l are specified, K and A? 7 n _j can 
be determined from equations (B3) and (B4). Generally, K will be a constant slightly 
greater than 1.0 in order to give smaller steps near the surface. The smaller steps near 
the surface increase the accuracy and efficiency of the computation for a typical turbulent 
boundary layer because the changes in mean velocity near the surface are usually much 
greater than those in the outer part of the boundary layer. Other parameters such as E 
and r may also change rapidly near the outer edge where, with K > 1.0, the step size 
could be considerably larger than that near the surface. Check solutions, however, with 
K = 1.00 and 1.02 have given essentially identical results for all quantities including E 
and Trji. Apparently, with K = 1.02, enough detail can be retained near the outer edge 
to avoid any loss of accuracy. 


Boundary and Initial Conditions 

The external velocity U e and its derivative dU e /dx^ must be specified functions 
of £. The nominal boundary -layer thickness 6, as used in equation (8), is taken as 
the point where 1 - F = which is a specified small constant (eg = 0.01, for example, 
in the computation of ref. 17). In order to insure an asymptotic boundary condition at the 
outer edge of the boundary layer, another boundary -layer thickness 6 e , defined as the 
point where 1 - F = e e , has been used with e e « eg. (For further discussion of edge 
boundary conditions, see section "Determination of Number of At] Steps.") For the 
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conventional boundary layers to be treated herein, both e e and eg would be small 

from equation (22) are 


\ (B5) 


profiles are specified at | Q (m = 1) 
from n = 1 to n = N, from which values of all variables are to be computed at the next 
station § 0 + or m = 2. If it is assumed that the initial profiles for E and V are 
locally similar (that is, 9V/d£, and 8E/a£ are neglected), these profiles may be com- 
puted from equations (27) and (29) and the input F 0 profile. In this procedure, equa- 
tion (29) is solved first for V and then equation (27) can be solved directly for E since 
M is a function of E from equation (19). 


positive numbers. The corresponding values of T] 

U c 


^6 


— 5 




v(2$ n 

n 


K24) 


The initial conditions in the form of F, E, and V 


Linear Finite -Difference Expressions 

The various types of derivatives in equations (27), (28), and (29) are replaced by 
linear difference quotients and equations (27) and (28) are evaluated at the intermediate 

grid points represented generally by the subscript m + i, n. The values at the interme- 

2 

diate points are computed to a satisfactory approximation as appropriate numerical aver- 
ages illustrated by the expressions: 


m + 


m + 


i ni 2 ( f m,n + f m+l,n,i-l) 

2 ’ ’ 

1 1 • 4 ( f m,n+l + *m,n + *m+l, n+l, i-1 + ^m+l,n,i-l) 

O’ r»> X 


2 2 


(B6) 


where f represents any desired quantity or function, and where the subscript i denotes 
the current step in the iteration cycle and the subscript i-1 denotes the value obtained 
in the previous iteration. The partial and total derivatives in equations (27) and (28) are 
then written with first-order accuracy as illustrated by the following examples: 



(B7) 


(B8) 
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8 i 

dF \ 

i 

M m+ 1 n+ l( F n»l."+l ' 
2’ 2 

_8t?I 

( 9 WJ 

1 A?7 + Ar? t 

m + ^r, n 11 11 x 

<1 


M 


m 


, 1 n _l( F m+l,n " F m+l,n-l) 


M 


2 2 


m + 


n + l( Fm ’ n+1 " F “»n) 


+ 


2-2 


M 


m 


A? ?n-1 

+ 1 n --( Fm » n ~ F m,n-l) 

O 5 O 


A7? n 


2’ 2 


At?. 


n-1 


(B9) 


All other terms and coefficients of derivatives appearing in the equations are 

linearized by evaluating them at the intermediate points m + ^ according to equa- 

* 2 ^ 

tions (B6). Furthermore, the (M - 1)( — ] term in equation (28) (which is treated as 

\ dr l I 

a linear equation for E) is evaluated as 


(“ ■ i) (if 


4 F _ i 


- F 


1 4 


= f M j - 1) 
1 \ m+-,n,i 

m+ ^,n,i \ ^ 


m + 2 ’ n + 2’ 1 m + 2’ n "? 1 


(AT? n + A ^) 5 


(BIO) 


An alternate expression, which is equivalent to the same accuracy as equations (B7) 
to (B9), is 


(M - 1) 


9F\ 2 " 

J 


= M 


A / r m+l,n+l,i-l J ' m+l,n-l,i-l \ / F m,n+1 ~ F m,n-A 


m + -,n,i 


m + -, n ,i 


A \ + A? Vl 


A7 ?n + A 7 n -1 / (Bll) 


In this manner the direct coupling between equations (27) and (28) is removed during any 
one iteration; however, the coupling is actually retained in the complete iteration proce- 
dure wherein the intermediate m + i values are continually improved, as indicated by 

u 

equations (B6), until convergence to the desired degree of accuracy is attained. 
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Computational Equations 

Mean momentum .- Substitution of the linear -difference quotients given by equa- 
tions (B7) to (B9) into equation (27) and collecting terms gives 

A n F m+l,n+l + B n F m+l,n + c n F m+l,n-l = D n ( B 


V 1 
m + 2> n 


m + I,n + | 


An 2(A77 n + A Vn _ x ) 2 A7j n A?j n _ 1 


3 = (2£) 2fi F + m + 2 >n+ 2 m + 2> n 

L J m + I m + ^ n 2 A7 7n A7 ?n-l 




C n — 


V 1 

m + ^n 


m + in-| 


2(Ar7 n + Arj^) 2 Arj n Arj n l 


1 _ . 1 


1 „ 1 

2 ’ 2 


D = _A F , , / f(2^) 2ff 1 m+ 2’ n+ 2 m+ 2’ n 

n n m,n+l ^ 2A„ n A Vl 


_c n F m n 1+ 

n m,n-l Ug j ^ 


m + 2 


Pa|] [1- F 1 2 

7 4 m 4 n 

m + 2 m +2 X 1 1 


1 + K 


k *- 2K 
1 + K 
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Since all quantities appearing in the A, B, C, and D coefficients are known 
from the completed calculation at the previous step m, or for the intermediate values at 

m + i from the previous iteration, equations (B12) represent a system of N - 2 equa- 

a 

tions with N - 2 unknowns which are the values of F m+ ^ from n = 2 to n = N - 1. 
^The boundary conditions require F m+1 ^ = 0 and F m+1 N = F e = 1 - e e ^. The matrix 
of the linear system (B12) is tri-diagonal; thus, the unknown F m+ j values are easily 
obtained by successive elimination of the unknowns (refs. 40 and 57) with the formulas 


■^131+ l,n ^n^m+l,n+l + *=n 


(B15) 


which are applied by starting the calculation at the outer edge of the boundary layer where 
Fm+1 n = B e an< ^ successively obtaining all other values of F m+ ^ down to the wall 
where F = 0 (eq. (20)). The functions G n and g n are computed from the recursion 
formulas 




G n = 


n 

B n + C n G n-l 


Sn 


D„ - C„g„ 
n n & n-l 

B„ + C„G„ i 
n nn-1 


(B 16 ) 


which are applied by starting at the wall (n = 1) and working out successively to 
n = N - 1. This procedure can be started by noting that with the wall boundary conditions, 
F m+ i j = 0, G-^ = g^ = 0 .from equation (B15). (See ref. 57.) Then, from equation (B16), 


G 


2 


A2 

B 2 





(B17) 


Determination of number of A 77 steps (N) .- Before the complete set of unknown 
values F m+ j can be obtained, it is necessary to determine the number of steps in 77 
required to satisfy the outer asymptotic boundary condition of F e = 1 - e e . The number 
of steps at any £ location is determined from the physical requirement that 



(B18) 
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where e g f is a small specified error criteria (generally a positive number for conven- 
tional boundary layers). The finite -difference form of inequality (B18) is written 


F m+1,N " F m+1,N-1 = Ar? N-l e e ( Bl9 ) 

Then from equation (B15) ^with F m+ ^ = F g ^ , 

F m+1,N-1 = G N-l F e + g N-l 

which upon substitution into inequality (B19) gives the requirement 


F e ( 1 _ G N-l) ' g N-l = A %-1 e e 


(B20) 


The computer program (appendix C) is written so that F e = F m _j N = 1 - e e ; that is, the 
value of e e is the same as that for the initial input profile. Thus with the error cri- 
teria e e and e e ' specified, successive values of G n and g n are used in inequal- 
ity (B20), and when the outer edge of the boundary layer is approached, the point at which 
the inequality is first satisfied determines the value of N - 1 and hence the number of 
A 77 steps required. 


If the value of r) e increases with | (with n constant), the computer program 
as described in appendix C is set up so that the required values of F m are obtained 

from the initial input profile of F as a function of 77 . Partly for this reason and partly 
because of the computing procedure used, it is necessary to specify F m _j N as a num- 
ber very close to unity (that is, e e ~ 1 X 10”^) in order to obtain ( — J = e e* as 

V 71 hr=ri e 


required by equations (B18) to (B20). 


Continuity .- The values of V ■. required in the A n and C n functions are 

m + |,n 

obtained from the continuity equation (29) which, evaluated in finite -difference form at 

point m + i n - - gives 
2 2 
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After a set of F m+ ^ n values are obtained, equation (B21) is used to compute an updated 

set of V . values where from the wall boundary conditions (eq. (20)) and the defini- 
m + -,n 

tion of V (eq. 26)), 


V 


m + |,i 




(B22) 


For zero mass transfer at the wall, iig = 0; otherwise, this normal velocity at the wall 
may be a specified function of x^. Successive application of equation (B21) from n = 1 
to n = N then gives the new values of V which may be utilized during the itera- 

m + 2’ n 

tion cycle for obtaining convergence of the F m+ ^ n values from equation (B12). In the 

present method, this convergence criteria is based on the shear stress from the 
expression 



5 e 


w 


(B23) 


where e w was generally specified as 0.01 unless otherwise stated. Since F w = 0, this 
inequality may be written in finite -difference form as 


F m+l,2,i 1 ~ F m+l,2,i-l ( e w) 


(B24) 


The convergence criteria e w is not an input quantity but appears in the program listing 
(appendix C) in the first statement in the section "Test for Iteration on FP." 

Since input values of U2 at £ Q (m = 1) are not generally available (for the com- 
putation of directly from eq. (26)), it is necessary to compute the V 

’ 1 + — ,n,i=l 

2 

values from equation (29) by dropping the £ derivatives. The resulting equation is 
written as 


V 


V 


1 + 2> n ’ i=1 


l.n 


= V 


l,n-l 


At vi r, 

2 L 


( 2 £) 2 n| 


^- l m=l 


( F l,n +F l,n-l) 


(B25) 
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For all the calculations presented, the value of n was retained as a constant. It 
was found that appreciable savings in computer time could be obtained by the use of a 
value of ii appropriate to the type of flow being computed. For example, for the flat- 
plate problem, if the value of ii at £ = 5 X 10 4 was changed from 0.5 to 0.8, the num- 
ber of A rj steps required at £ = 10® was reduced from 290 to 125. Actually, the 
problem had to be reinitialized at | = 5 x 10 4 to accommodate the new value of ii, but 
the computer time from that point on was reduced by approximately one-half. 

Energy .- The finite -difference form of equation (28) is derived in the same way as 
equation (B12) and, with the use of equation (Bll), is written in the same form as equa- 
tion (B12) 

A n E m+l,n+l + E n E m+l,n + ^n E m+l,n-l = E n (B26) 

AAA A 

The A n , B n , C n , and D n coefficients are obtained as 


v 1 
m+ 2’ n 


2 (A 7 , n + ATJn-i) 


KD j i 

m + 2’ n+ 2 
2 A rj n Ari n _i 


B n = 


(2 Q 


2n 


A? 


KD 


1 + K+E L . 1 


F 1 + 
m + I m+ 2’ n 

Li 


m + _ n + - m + -,n-- 

2 At j n A77 n-1 


m + 2 ,n 


K*D 1 1 

m+ 2’ n -2 


2 (AT? n + Aij^j) 2AT 7n AT ?n . 1 


D n — A n E m,n+l + 


" C n E m,n-l ' 2 


(Ml 


2n 


AS 


1M)1 


KD , . + K*D . . 

m+^n+g- m + i n - 2 


1 m + i,n 


2 A7 ?n Ar hi-l 


E r 


U e A? 


dU c 


■ A« 


E F 

1 _ „ . 1 


l\d? / i m + An m + An 
m + 2 m + 2 2’ 2 ’ 


(B27) 


M 1-1 
m + £,n 


(AT7 n + At 7ii _ 1 ) : 


| , ( E m,n+l " F m,n-l) ( F m+l,n+l F m+l,n-l) i - C 


D i E i 
m + ^,n m+g.n 

i 2 <Pn 2 

6,m + - 
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It should be noted that in these computations, 

= •<{%) 

because of the definition of 77 (eq. (22)), The procedure for solving the system of equa- 
tion (B26) is identical to that outlined, except that the number of Atj steps required is 
not recomputed separately to satisfy directly the edge physical requirement that 

| approach 0. Instead, the value of N already determined from equation (B20) is 
e 

used also for the solution of the energy equation, and at this value of N, the edge boundary 
condition of Ej^ = E e is imposed where E e is generally specified as a small finite 

number. It has been found that in all cases computed, this procedure gives 
approaching 0 to the desired degree of accuracy. 



/9E 
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DIGITAL COMPUTER PROGRAM DESCRIPTION 

By Carolyn C. Thomas 
Langley Research Center 

This program computes the nonsimilar development of an incompressible turbulent 
or laminar boundary layer by the finite -difference procedure described in appendix B. 

Minimum machine requirements on the Control Data 6400/6600 computer system 
are 67000 octal locations of core storage. The program is written in FORTRAN 2.0 
which is compatible with FORTRAN IV in most instances. The time required to calcu- 
late a grid point is approximately 0.0014 second per iteration. The total computing time 
required for a typical Goldberg case was 122.5 seconds where the pertinent inputs were: 
K = 1.02, Arjj = 0.005, n = 0.5, e w = 0.01, e e ’ = 1 x 10" 5 , and A£ = 1 x 10 3 . The 
average number of A 77 steps across the boundary layer was about 95 and there were a 
total of 910 A£ steps for this particular case. With the iteration criteria of % = 0.01, 
the solution was iterated only 34 times over 910 ^-stations (or columns) computed. Sev- 
eral of these iterations were required in the first 3 or 4 A£ steps where the initial 
input profiles for F and E are being adjusted to achieve consistency with the equa- 
tions of motion. 


Input 

One card of case identification using all 80 columns. The remaining input is loaded 
by using the FORTRAN version 2.0 NAMELIST. The input symbols are as follows: 


Symbol (see section "SYMBOLS") Machine name Comments 

f FPTAB Input table, of up to 600 values 


F = F(r?) 


NUMFP Number of values in FPTAB table 

(Integer value) 

ETATABF ETA table that corresponds to FPTAB 


Eo 

P 


EO 

RHO 

NU 
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a 

ALPHA 


K 

KN 

Ratio of Ar) at n + 1 to Arj at 

R 0 

RO 


C 

C 


K 

K 


£ 0 

XIO 


n 

NBAR 


No 

NMAX 

Initial number of steps in boundary 
layer (Integer value) 

A4 

DELXI 


At) 

DELETA 

Initial value of A r q 

v*/v e 

ETAMOE 


x 2 /5 

TX20DEL?\ 

Table of up to 8 values 

1/5 

TLODEL J 


e e 

EPSLON 



NUMETA 

Total number of values expected at 
end of case (can be 1 to 600) 
(Integer value) 


MPRINT 

Table of values at which complete 
print out is desired - up to 
40 values 


FINLPRT 

Final value from MPRINT table 
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£ 


U e 


dUe/d£ 


E = E(rj) 


X3TAB Table of values used to obtain U e 

value from UETAB and dU e /d£ 
value from DUETAB (50 values) 


UETAB Table of up to 50 values. If U e is 

constant, all values in table should 
be set to this value. 


DUETAB 


(ETAB 


INUME 


E TAT ABE 


Table of dUe/d£ values (up to 
50 values) 

Input table of up to 600 values 

Number of values in ETAB (Integer 
value) 

ETA table that corresponds to ETAB 
input 


IEREAD =0 (E table read), ^0 (E table com- 
puted) (Integer value) 


The program will interpolate between input values of F and E to fill a table having 
NUMETA values with a spacing of Ary. As mentioned in appendix B, it is necessary to 
specify Nq at a point on the input profile where 1 - F n < 1 x 10 in order to obtain 
the desired slope Hp-) = e e ' at subsequent profiles. Furthermore, because of the 
V v /ri=v e 

interpolation procedure used to fill the F table, the last value of F in the input table 
must be slightly greater than F at n = N. 

Output 

Identifications of case and input data are printed at beginning of output. 


Symbol 



Machine name 


Comments 


MCOL (designates m column)' 


XIA 


These are printed at 
every m column 
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s 

XI 

*\ 


Cf 

CF 



9 

THETA (ft) 





> 

5* 

DELSTAR (ft) 


H* 

H 



R 6 

RTHETA 

> 


V 

ETA 



F 

FP (computed column) 

E 

E (computed column) 

V 1 
m+ 2’ n 

VA 



Uj/ll* 

UUSTAR 



u*yA 

USYNU 



U e -U! 

UUEUS 



u* 




n 

KCH 


J 

/ _\2 

> 



< M - «(lf) 

PROD 



CDE/tjS 2 ^ 2 

DISIP 

> 


(2£) 2fi F + 

CONVEC 



JLl d 
917 \ 917/ 

DIFUS ^ 




These are printed at 
every m column 


Group I 


Group II 
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2 ( 2 £) 


2n du e FE 
d£ U e 


r^j2pe 


DUETRM 

TORHOE 


Group II 


Print out is controlled by MPRINT table. Groups I and II are printed only when f 
corresponds to a value in MPRINT table. XI table must be filled with values in ascending 
order. 
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NASA -LANG LEY RESEARCH CENTER 


02 2b 

CATEGORY 

02 27 

LANGUAGE 
NO 1 

02 L2 

LANGUAGE 
NO 2 

fin 

Ff5R6M 



01 4 

LAR 


01 7 PROGRAM NO 

D2170 


COMPUTER PROGRAM ABSTRACT 


01 20 TITLE OF PROGRAM f61 CHARACTERS MAX«VUV 

Calculation of Turbulent and/or Laminar Boundary Layers 


01 14 DATE 

072268 

PARENT PROGRAM 


02 18 PROGRAM NO 


02 3? KEY WORDS 8 MAXIMUM SEPARATED BY COMMAS' 


Bound ary Layer, Turbulent, Laminar, Finite Difference 


Proce dure, Implicit 


WHO TO CONTACT ABOUT THE PROGRAM 


05 14 CONTACT 

C. C. Thomas 


05 50 INlTlATEO 


0966 


DATES 

f05 54 COMPLF Tf D 


CARO NUMBER 


05 28 SITt 

LAR 


05 31 ORGN CODE 

11.160 


P5 35 PROJECT NO 

RGL 151 


05 45 
NASA 
CENTER 


_ .lb St REVISION CODE 

I \ A REVISION 

B CANCELLATION 


05 50 MANMONTHS 


5« 60 61 62 63 


05 48 STATUS 

CD A UNDER DEVELOPMENT 
B OPERATIONAL 
1 I C COMPLETED 

TIME AND COST FOR DEVELOPMENT 
I MACHINE 1 05 65 COMPUTER TYPE 


6000 


ELITE MARGIN 


ABSTRACT 


An implicit finite -difference scheme is used to compute 
incompressible boundary layers with arbitrary pressure 
gradient and wall blowing. To compute a turbulent boundary 
layer, an eddy viscosity concept is used wherein the eddy 
viscosity is assumed to be a certain function of the kinetic 
energy of turbulent fluctuations and the scale of the turbu- 
lence. The differential equation for this turbulent energy is 
then solved simultaneously with the conventional momentum 
equation to obtain the distribution of the mean velocity and the 
distribution of the turbulent energy. Auxiliary functions are 
supplied for the eddy viscosity and the turbulence scale. The 
objectives of the present investigation are to develop the com- 
puting procedure and to determine the range of applicability of 
the auxiliary functions. 


□ A THIS PROGRAM j 
IS NOT EOR 
SHARING 


05 74 TOTAL COST 
■DOLLARS 
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24 

26 

26 

27 

28 

29 
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Flow Diagram 


Start \ 



Read and Print 
input data 

Initialization 


Compute ETA 
table — Compute 
FP and E tables 
by using FPTAB 
and ETAB tables 
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Program Listing 


JOB.01 ,2000*67000. D2 1 70 , 3 1 366 * 1 , CAROL YN C THOMAS , , 1 1 92-C Rl6i 

L I NECNT (40000) 

RUN (S ) 

SETINDF. 

LGO • 

EXIT. 

DMP < FL ) 

PROGRAM D2170 ( INPUT , OUTPUT ,TAPE5=INPUT , TAPE6 = OUTPUT ) 

********** 

PROGRAM TO COMPUTE THE NON-SIMILAR DEVELOPMENT OF AN INCOMPRESSIBLE 
TURBULENT OR LAMINAR BOUNDARY LAYER 
********** 

*****CAROLYN C. THOMAS*****ANALYSI S AND COMPUTATION D I V I S I ON***** 1 967 

D I MENS ION ID (8 ) .FP (600 . 2 ) , FP1 ( 60 0 ) ,FP2 < 600 ) , FPTAB (600 > »FPA (600 ) . 

1 FPA 1 (600 ) iE (600.2 ) . E A ( 600 ) , ET AB (600 ) iETA (600 ) .ETATABF ( 600 ) , 

2 ETATABE <600 ) ,M <600 »2 ) .MAI (600 ) »MA2 <600 ) »MA <600 ) , D<600,2 ) , DA (600 ) » 

3 DA 1 <600 ) , DA2 <600).A(600),B(600>*T(600)»S<600),Q(600),U(600), 

4 W <600 ) iZ (60C ) , VA ( 600 ) , VAI NT <600 ) . X2 (600 ) . TX20DEL (8 ) , TLODEL (8 ) . 

5 LODEL (600.2). LODELA ( 600 ) . DELT A (2 ) . CAPG (600 ) . SMLG ( 600 ) . UET AB ( 50 ) . 

6 DUETAB ( 50 ) . XI TAB (50 ) . MPRI NT (40). DEL ETA (600 ) 

EQUIVALENCE ( M A 2 < 1 ) .MAI (2 ) ) . (DA2 ( 1 ) . DAI (2 ) ) 

EQUIVALENCE (FP ! ( 1 ) ,FP ( 1 i 1 ) ) i ( FP2 ( 1 ) ,FP( 1 .2) ) 

EQUIVALENCE (A(l).Q(l)).(d(l),U(l)).(T(l).W(l)).(S(I).Z(l)) 
EQUIVALENCE (FPTAB ( 1 ), FP ( 1 . 2 )).( ETAB ( 1 ), E ( 1 . 2 )> . 

1 ( ETATABF ( 1 ) . VA I NT ( 1 ) ) . (ETATABE ( 1 > .M ( 1 ,2 > ) 

REAL L. LODEL, M,K. MAI , MA2 , NU , NBAR , MA , LODELA , MPR I NT , IFPT.KN 
LOGICAL KFIRST 

NAMELIST /DA T A/FPT AB, ET AB, ETAT ABF, ET ATABE, EO.RHO.NU, ALPHA, RO.C.K, 

1 X I 0,NBAR,NMAX , DELX I , DELET A , ET AMOE , EE , T X20DEL , TLODEL , EPSLON , 

? NUMETA , MPRI NT ,F I NLPRT , X I T AB , UETAB , DUET AB , I ERE AD , NUMFP , NUME , KN 

INPUT AND INITIALIZATION 

NUMF=0 

1 READ (5,2) ID 

2 FORMAT ( 8 A 1 0 ) 

READ ( 5 » DAT A ) 

I FPT=MPR I NT ( 1 ) +DELX I 
IP=1 

WRITE (6, DATA) 

PRINT 3, ID 

3 FORMAT ( 1 HI //8A1 0 /// ) 

UEPRES = UET AB ( 1 ) 

DELUEPS = DUETAB ( 1 ) *DELX I 
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MCOL= 1 
I TER = 0 
XI =XIO 

X I A= (XI O+DELXI +XI 0 )/2. 0 
J= 1 

KF I RST= • TRUE • 

XMULT1 =2.0*KN/ ( 1 .O + KN ) 

XMULT2=2.0/< 1 .O+KN) 

NEM1 =NUMETA-1 

TW0NBAR=2. 0*NBAR 

NMAX20=NMAX+20 

NMAX40=NMAX+40 

MORDF = 2 

MORDE=2 

ETA (1 ) =0 .0 

IF (IEREAD .NE. 0) GO TO 200 
NECK=.90*FLOAT (NUMF ) 

ETANECK = ETATABE (NECK ) 

DO 300 ITAB=1 iNUMETA 

CALL FTLUP ( ETA ( I TAB )iFP(ITABil ) » MORDF . NUMFPiETATABF ( 1 ) « FPTAB ( 1 ) ) 
IF (FP(ITAB.I) .GE. .9) MORDF = 1 
IF (IEREAD .NE. 0) GO TO 250 

CALL FTLUP (ETA ( I TAB ) i E ( I TABt 1 ) » MORDE » NUME i ETATABE ( 1 ) » ETAB ( 1 )) 

IF (ETA(ITAB) . GE . ETANECK ) M0RDE= 1 
250 DELETA ( I TAB+1 ) =KN*DELETA ( I TAB ) 

ETA ( I TAB+1 ) = ET A ( I TAB )+ DELETA < I TAB > 

300 CONTINUE 

IF (IEREAD .EQ. 0) GO TO 400 

COMPUTE E(N.l) FOR N = 1 i NUMET A (INITIAL VALUES) 

ET AE = ET A (NMAX ) 

ETAM=ETAMOE*ETAE 
DO 4 N = 1 » NUMET A 
ETAQ=ETA (Nl/ETAM 

4 E (N. 1 ) =EO*ETAQ**2* ( EXP ( . 5* ( 1 . 0-ETAQ**2 ) ) >**2 
400 CONTINUE 
VA ( 1 ) =0.0 
DO 126 NA= 1 ♦ NEM 1 
VAO= VA ( NA ) 

VA (NA+1 ) =V AO— DELETA (NA ) *TWONBAR* ( 2.*X I A >** ( TWONBAR- 1.0) *FP (NA + 1.1 ) 
126 CONTINUE 

FPTEST = FP ( NMAX , 1 ) 

COMPUTATION OF DELTA AND DELTA (E) 
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39 XFP=.990 
393 NEWE=0 

IF (J .NE. 1 ) GO TO 391 

CALL FTLUP < XFP', XET A , I , NMAX » FP 1 ( 1 ) » ET A ( 1 ) ) 

GO TO 392 

391 CALL FTLUP ( XFP tXETAi 1 . NMAX » FP2 ( 1 ) » ET A ( 1 ) ) 

392 DELTA ( J ) = ( NU* ( 2 . 0*X I > **NBAR ) /uEPRES*XETA 

C FREQUENTLY USED EXPRESSIONS. 

NMAX20=NMAX+20 

IF (NMAX20 .GT. NEM 1 > NMAX20=NEM1 
NMA'X40 = NMAX+40 

IF ( NMAX40 .GT. NEM 1 ) NMAX40=NEM1 
TWOXIA=2.0*XIA 
TX I ANBR = TWOX I A**NBAR 
TX I A2NB=TW0X I A**TWONBAR 
C 

C COMPUTATION OF AUXILIARY FUNCTIONS R,H(R)*M,DiL 

C 

HR=HKR=0 • 0 
NMAX71 =NMAX40+3 1 

IF (NMAX71 .GT.NUMETA) NMAX71 =NUMETA 

40 DO 6 NROW=l .NMAX71 

X2 (NROW ) = (NU*TX I ANBR ) /UEPRES*ETA (NROW ) 

X20DEL=X2 (NROW ) /DELTA ( J ) 

IF ( X20DEL .LE. 1.4) GO TO 46 
LODEL ( NROW » J >=0.01 
GO TO 47 

46 CALL FTLUP ( X20DEL , LODEL ( NROW . J >• 1 . 8 . TX20DEL ( 1 ). TLODEL ( 1 ) ) 

47 L = LODEL (NROW , J ) *DELTA ( J ) 

IF (E(NROW.J) .LT. 0.0) E ( NROW « J > =0 . 0 
R=SQRT ( E (NROW, J ) ) *UEPRES*L/NU 
RRO=R/RO 

IF (HR.EQ.l . )GO TO 45 
I F ( RRO .LT. 1.25) GO TO 41 
HR= 1 .0 
GO TO 45 

41 IF (RRO .LT. 0.75) GO TO 42 
HR=RRO- (RRO- 0.75 )**2 

GO TO 45 

42 IF (RRO .GE. 0.0) GO TO 44 
PRINT 4 3 , R , RO , NRO W 

43 FORMAT (//32H M=l,R/RO NEGAT I VE , WHAT HAPPENED/3H R=E16.8,5X, 
1 3HR0 = E16.8 ,5X, 2HN= I 5/ ) 

44 HR =RRO 
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45 EPSLONR=ALPHA*R*HR 

M ( NROW i J)=l . O+EPSLONR 

RK=R*K 

RKRO=RK/RO 

IF (HKR.EQ. 1 . )GO TO 50 
IF ( RKRO ,LT. 1.25) GO TO 451 
HKR= 1 .0 
GO TO 50 

451 IF (RKRO .LT. 0.75) GO TO 452 
HKR=RKRO- (RKRO-O.75 )**2 

GO TO 50 

452 IF (RKRO .GE. 0.0) GO TO 49 
PRINT 48 i RK iROi NROW 

48 FORMAT (//33H M=1,RK/R0 NEGAT I VEi WHAT HAPPENED/4H RK=E16.8.5X« 
I 3HRO=El 6« 8 t5X i 2HN= 15/ ) 

49 HKR =RKRO 

50 EPSLNKR=ALPHA*R*K*HKR 
D ( NROW . J ) = 1 . O + EPSLNKR 

6 CONTINUE 

X I BAR = TX I A2NB/DELX I 
IF (KF I RST ) 61.51 

SET INITIAL VALUES OF VARIABLES IN COLUMN 2 EQUAL TO VALUES IN 

61 DELTA (2 )=DELTA ( 1 ) 

DO 5 N= 1 .NMAX71 
LODEL (N.2) =LODEL (N, 1 ) 

M (N ,2 ) =M (N, 1 ) 

5 D ( N.2 ) =D (N. 1 ) 

DO 53 N = 1 .NUMETA 
FP (N, 2 ) =FP (N . 1 ) 

E ( N . 2 )=E (N. 1 ) 

VA 1 NT ( N ) — VA ( N ) 

53 CONTINUE 

INITIALIZE X I . UE . AND DELUE FOR EACH COLUMN 
SLOPE = FP (2.2 J/DELETA ( 1 ) 

J = 2 

X I =X I +DELX I 

TWOX I =2 • 0*X I 

TXINBAR=TWOXI**NBAR 

UEPREV=UEPRES 

DELUEPV=DELUEPS 

CALL FTLUP (XI. UEPRES . 1 .50. XITAB ( 1 ) . UET AB ( 1 ) ) 


COLUMN 1 
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CALL FTLUP (XI i DELUEPS ill 50 iXIT Ab ( 1 > * DUET AB ( 1 ) > 

DELUEPS=DELUEPS*DELX I 

COMPUTATION OF AVERAGE VALUES OF M.FP.UE.AND DELUE FOR MOMENTUM EQUATION 


51 FPA 1 =0 • 0 

MAI ( 1 ) = 1 .0 

FPA ( 1 )= <FP ( 1 . 1 )+FP( 1 ,2 ) 1/2.0 
DO 7 NAM . NMAX40 

FPA (NA+I 1 = (FP (NA+1 » 1 1 +FP( NA+1 , 2 1 1/2.0 

FPA1 (NA+1 )= ( FP (NA , 1 ) +FP (NAi 2 I+Fp (NA + 1 ,1 ) +FP (NA + 1 ,2 ) ) /4.0 
MA 1 (NA+1 ) = (M (NA , 1 ) +M ( NA * 2 ) +M ( NA+ 1.1 )+M ( NA+1 ,2 ) )/4.0 

7 CONTINUE 

UEA= (UEPRES+UEPREV )/2. 0 
DELUEA= (DELUEPS+DELUFPV 1/2.0 
UEA1 =UEA 
DELUEA1 =DELUEA 

KF IRST= .FALSE. 

COMPUTATIONS TO OBTAIN G CONSTANTS FOR MOMENTUM EQUATION 

NMAX39=NMAX40- 1 
76 DO 8 KG=2,NMAX39 

CON7=2.0*DELETA ( KG ) *DELETA ( KG- 1 ) 

CONI =VA (KG )/ (2.0* (DtLETA (KG l+DELETA (KG-1 ) > ) 

C0N2 = (MAI (KG )*XMULT1+MA2(KG)*XMULT2 1/C0N7 
C0N3=X IBAR*FPA (KG ) 

A (KG )=C0N1 -MA2 (KG ) *XMUL T2/C0N7 
B ( KG )=CON3+CON2 

T (KG ) =- (CONI +MA 1 ( KG ) *XMULT 1 /C0N7 ) 

C0N4=C0N3— C0N2 

C0N5=X IBAR+DELUEA/UEA* ( 1 .0-FPA (KG )**2 ) 

S (KG ) =-A (KG ) *FP (KG+1 , 1 ) +CON4*FP (KG. 1 )-T (KG ) *FP (KG-1 . 1 J+CON5 

8 CONTINUE 

CAPG ( 2 ) = — A (2 >/B (2 ) 

SMLG (2 )= (S (2 )~T (2 )*FP (1 ,2 ) )/B (2 ) 

DO 9 KG=3,NMAX39 

CAPG (KG ) =-A (KG )/ ( B (KG >+T (KG ) *CAPG (KG-1 ) > 

9 SMLG (KG ) = (S (KG ) -T (KG ) *SMLG (KG-1 ) ) / (B (KG ) + T (KG')*CAPG ( KG-1 ) ) 

COMPUTE NEW BOUNDARY HEIGHT 
KMAX=NMAX- 1 0 
DO 90 JJ=KMAX. NMAX20 

TSTVAL=EPSLON*DELETA ( JJ ) 
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CKVAL=FPTEST* ( 1 .O-CAPG ( J J ) ) -SMLGUJ > 

NMAX= JJ+2 

IF ( CKVAL .LE. TSTVAL ) GO TO 91 

90 CONTINUE 

COMPUTE NEW FP AND VA VALUES 

91 NRACK=— (NMAX-1 ) 

FPPREV=FP (2 ,2 ) 

DO 10 NF = NB ACK , — 2 
KF= I ABS (NF ) 

FP (KFi2 ) =CAPG(KF )*FP (KF + 1 i2 J+SMLG (KF ) 

10 CONTINUE 
C 

VA ( 1 )=0.0 
DO 18 NA = 1 * NMAX40 
52 VAO = VA ( NA ) 

18 VAINA + l ) =VAO— DELETA (NA ) /2 • 0*X I BAR* (FP (NA+ 1,2) -FP (NA+ 1,1) +FP (NA , 2 ) - 
1 FP ( NA , 1 ) ) -DELETA ( NA ) *TWONBAR*T WOX I A** ( TWONBAR— 1 . 0 ) *FPA 1 (NA+ 1 ) 

C 

C COMPUTE CF AND THETA 

1003 SLOPRE V=SLOPE 

SL0PE=FP (2,2 ) /DELETA ( 1 ) 

CF =2. 0/TX INBAR* SLOPE 
PRINT 1 01 2 , SLOPREV, SLOPE 

1012 FORMAT </9H SLOPREV= E 1 6 .8 , 5X , 6HSL0PE =E 1 6 . 8/ ) 

DFLST AR = 0 . 0 
THETA=0.0 

DO 1020 KGT = 1 , NMAX 

DELSTAR = DELSTAR+ ( 1 ,0-FP (KGT ,2 ) ) *DELET A (KGT ) 

10 20 THETA = THFTA + FP (KGT,2)*(1.0-FP(KGT,2) ) *DELET A ( KGT ) 

DELSTAR=TX 1 NB A H*NU/UF A* DELS T A R 
THETA=NU*TX I NB f R/UE A*THETA 
RT HET A = UE A *T HET A/NU 
H=DELSTAR/THETA 

PRINT 1 050 , MCOL , X 1 A , X I , CF , T HET A , DELS T A R , H , RT HE T A 
1050 FORMAT ( 1 9H READY FOR COLUMN I 5/5H X I A = E 1 6 , 8 , 5X , 3 HX I = E 1 6 .'8 , 5X , 

1 3HCF = E 1 6.8, 5X , 6HTHETA = E1 6.8/ 1 OX, 1 1 HDELTA ST AR = E 1 6 . 8 , 5X , 

2 2HH = E16.8,5X,7HRTHETA = E16.8/ ) 

C 

C TEST FOR ITERATION ON FP 

I F ( ABS ( ( FP ( 2 , 2 ) -FPPREV ) /FPPREV ) -. 005 . GT . 0.0) GO TO 1005 
C 

IF ( ABS (X I-IFPT ) .GE. DELXI) GO TO 800 
I P= I P+1 
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I FPT=MPR I NT ( IP ) 

92 WRITE (6,999) 

999 FORMAT ( 9X , 8HET A ( KCH ) , 8X , 9HFP ( KCH , 2 ) , 9X , BHE ( KCH , 2 ) , 1 OX , 

1 7HVA (KCH ) ,6X, 1 1HUUSTAR (KCH ) *7X, 1 OHUSYNU (KCH > , 7X , 1 OHUUlUS ( KCH ) , 

2 5X.3HKCH) 

SQRTCF2 = SQRT (CF/2.0 ) 

DO 1002 KCH=1_,NMAX 
UUSTAR=FP (KCH, 2 )/SQRTCF2 
USYNU=ETA (KCH ) *TX INBAR* SORT CF2 
UUE(JS= ( 1 ,0-FP(KCH,2 ) )/SQRTCF2 

PRINT 1 001 ,ETA (KCH ) ,FP ( KCH, 2 ) , E ( KC H , 2 ) , V A ( KCH ) , UUST AR , US YNU , 

1 UUEUS , KCH 

1001 FORMAT (7E17,6, 18) 

1002 CONTINUE 

COMPUTE AND PRINT PROD , D I S I P , C ONVEC , D 1 FUS , DUE TERM , AND TORHOE 
DO 16 KGE = 2 , NM A X 

CON7 = 2. 0*DELETA (KGE ) *DELET A ( KGE — 1 ) 

F0RDETA = 2 .0* (DELETA (KGE ) +DELET A (KGE— 1 ) ) 

CON 1 2 = EA (KGE J/LODEL A (KGE )**2’ 

PROD= ( MA (KGE )-!.)/( DELETA (KGE )+ DELETA ( KGE— 1 ) ) **2* (FP (KGE + 1 ,1 )- 
1 FP (KGE- 1,1 ) ) * (FP (KGE+1 , 2 ) -FP (KGE- 1 , 2 ) ) 

D I S IP=TX I A2NB*C* (NU/ ( UEA*DELTAA ) )**2*DA (KGE )*C0N1 2 

CONVEC=VA (KGE ) /FORDETA* (E (KGE+1 ,2 )+E ( KGE+ 1 , 1 ) )-VA ( KGE ) /FORDETA* 

1 ( E (KGE- 1 ,2 )+E ( KGE- 1 , 1 ) )+TXI A2NB/DELX I *FPA ( KGE ) * ( E ( KGE , 2 ) -E (KGE, 1 ) ) 
D I FUS = DA2 (KGE ) *XMUL T2/C0N7* ( E ( KGE+ 1 ,2 )-E (KGE , 2 )-E (KGE * 1 ) + 

1 E ( KGE+ 1,1 ) )+DA 1 (KGE ) *XMULT 1 /C0N7* (E (KGE- 1 ,2 ) -E (KGE , 2 ) -E (KGE, 1 ) + 

2 E (KGE- 1,1')) 

DUETRM=2,0*TXI A2NB/UEA+DELUEPS/DELX I *EA (KGE )*FPA (KGE ) 

IF (EA(KGE) .NE. 0.0) GO TO 1008 
T ORHOF = 0.0 
GO TO 1009 

1 0 08 TORHOE = (MA(KGE)-l .0)/(EA(KGE)*2.*TXIANBR)*l . /FORDETA* 

1 ( (FP (KGE+1 ,2 ) +FP (KGE+ 1 , 1 ) )- (FP (KGE-1 , 2 ) +FP ( KGE— 1 ,1 ) ) ) 

1009 IF (KGE .NE. 2) GO TO 705 
PRINT 706 

706 FORMAT (/I OX , 1 OHPRODUCT I ON , 9X , 1 1 HD 1 SS I PAT I ON , 1 OX , 1 OHCONVECT I ON, 

1 1 1 X,9HD IFUS ION, 1 2X,8HDUE TERM , 1 2X , 8HTAU0RH0E ) 

705 PRINT 707, PROD, D IS I P , CONVEC * D I FUS , DUETRM , TORHOE, KGE 

707 FORMAT (6E20.8,I8) 

16 CONTINUE 

800 MC0L=MC0L+1 
NFWE = 1 
GO TO 105 
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C INITIALIZATION FOR NEXT COLUMN 

103 DO 102 KS=1 i NMAX20 
FP (KS . 1 ) =FP (KS .2 ) 

102 E (KS, 1 ) =E (KS ,2 ) 

106 X I A = X I A + DFLX I 
KFIRST=.TRUE. 

J= 1 

I TER=0 

IF (XI .GE. FIN!. "RT) GO TO 1 
GO TO 39 

CHECK TO SEE IF MAXIMUM ITERATIONS HAVE TAKEN PLACE 
1005 I TER= I TER + 1 

IF (ITER .LT. 10) GO TO 105 

HALVE THE INTERVAL DELX I AND RE- I N I T I AL I ZE AND TRY AGAIN 
I TER=0 
J= 1 

XI =X I —DELX I 
DELXI=DELXl/2.0 
X I A= ( X I+X I +DELX I ) /2.0 
KF I RST = • TRUE. 

UFPRFS=UFPRE V 
DELUEPS=DELUEPV 
DO 104 NVA=1 ,NMAX40 

104 VA (NVA ) =VA INT (NVA ) 

GO TO 39 

COMPUTATION OF AVERAGE VALUES OF EiMi AND D FOR ENFRGY EQUATION 

1 05 DELTA A= (DELTA (1 ) +DELT A ( 2 ) )/2.0 
DAI (1 ) = 1 .0 

DO 11 NA = 1 « NMAX40 

DA (NA ) = (D (NA , 1 ) +D ( N A , 2 ) ) /2 . 0 

DAI (NA+1 ) = (D (NA i 1 ) +D < NA . 2 > +D ( NA+ 1 » 1 )+D(NA+l , 2 ) ) /4 . 0 
DAI (NA+1 ) = DA 1 (NA+1 )*3.0 

LODELA (NA ) = (LODEL (NA, 1 J+LODEL (NA,2 ) ) /2 . 0 
EA (NA ) = (E (NA , 1 ) +E ( N A , 2 ) ) /2 . 0 
11 MA (NA > = (M (NA . 1 ) +M (NA, 2 ) ) /2 .0 

COMPUTATIONS TO OBTAIN G CONSTANTS FOR ENERGY EQUATION 

DO 12 KGE = 2 , NM A X 

CON7=2.0*DELETA (KGE )*DELETA ( KGE— 1 ) 

CON6 = VA (KGE )/ (2.0* ( DELE T A ( K GF ) +DELETA (KGE— 1 ) ) ) 
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Q < KGF ) =C0N6-DA2 (KGE ) *XMULT2/C0N7 
C0N8=X I BAR*FPA ( KGE ) 

C0N9= (DAI (KGE )*XMULT1 +DA2 (KGE )*XMULT2 )/C0N7 
U (KGE ) =C0N8+C0N9 

W (KGE >=- (CON6+DA1 (KGE ) *XMULT 1 /C0N7 ) 

IF (LODELA'(KGE > »NE. 0.0) GO TO 113 
CONI 2 = 0.0 
GO TO 12 

113 CONI 2=EA (KGE )/LODELA ( KGE ) **2 

Z (KGE ) =-Q (KGE )*E ( KGE+ 1 , 1 ) + ( C0N8-C0N9 )*E (KGE i 1 ) - W (KGE ) *E (KGE-1 . 1 ) — 
12.0*XIBAR*DELUEA/UEA*EA (KGE)*FPA (KGE ) -X I 8AR#DELX I *C* ( NU/ ( UEA * 
2DELT A A ) ) **2*DA ( KGE >#C0N12+ ( (MA (KGE ) - 1 .0 )/ (DELETA (KGE ) + 

3LELETA (KGE-1 ) )**2 

3 * (FP (KGE + 1 . 1 ) — FP ( KGE— 1 . 1 ) ) * (FP ( KGE+ 1 i 2 ) -FP (KGE— 1 ,2 ) ) ) 

12 CONTINUE 

CAPG (2 ) = — Q (2 )/U (2 ) 

SMLG ( 2 ) = ( Z ( 2 ) — W ( 2 ) *E ( 1 » 2 ) )/U(2) 

DO 13 KGE = 3 « NMAX 

CAPG (KGE ) =-Q (KGE ) / (U ( KGE >+W (KGE )*CAPG (KGE-1 ) ) 

1 3 SMLG (KGE ) = (Z (KGE ) -W (KGE )*SMLG (KGE-1 ) ) / ( U ( KGE )+W (KGE ) *CAPG (KGE-1 ) ) 
C 

C COMPUTATION OF NEW E VALUES 

C 

DO 14 NE=NBACK . — 2 
KE= I ABS ( NE ) 

14 E (KE.2 )=CAPG (KE ) *E (KF+1 «2>+SMLG(KE) 


C 

C 

C 

1922 


CHECK TO 
IF 

IF (NEWE 
GO TO 39 
END 


SEE IF ITERATING OR IF 
ITERATING GO TO 39 *** 
.EQ . 1 ) GO TO 103 


YOU HAVE ACCEPTABLE ANSWERS 
IF ACCEPTABLE ANSWERS GO TO 


103 
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TABLE I.- FUNCTIONS FOR MEAN SCALES 
OF TURBULENCE 


x 2 /6 


1/6 


^ 0.33 

^ 0.20 

^ 0.25 

0 

0 

0 

0 

.2 

.20 

.20 

.20 

.4 

.30 

.20 

.25 

.5 

.33 

.20 

.25 

.6 

.32 

.20 

.25 

.7 

.30 

.20 

.25 

.8 

.26 

.20 

.20 

= 1.4 

.01 

.01 

.01 



TABLE EL- INITIAL VELOCITY AND ENERGY PROFILES AND OTHER 
INPUT DATA FOR ADVERSE PRESSURE GRADIENT CASE 


[Values based on Goldberg's data (ref. 43): £ 0 = 8 X 10^; 

v = 1.8 X 10“ 4 ft2/sec (16.7 x 10-6 m 2/ s ) ; p = 0.0765 lb/ft 3 
(1.23 kg/m 3 ); U e>0 = 77.8 ft/sec (23.7 m/s); A£ = 625; 

= 0.025; n = 0.5; e ej0 = 9 x 10 " 4 5 e e ' = 1 X 10“ 5 ; 
e w = 0.01; ^=lx 10-2] 


V 

F ?o 

Et 

so 

0 

0 

0 

.1 

.215 

.22 X 10" 2 

.2 

.415 

.54 

.3 

.540 

1.05 

.4 

.580 

1.04 

.5 

.608 

.92 

.6 

.626 

.86 

.8 

.644 

.78 

1.0 

.662 

.72 

1.2 

.681 

.675 

1.4 

.700 

.64 

1.8 

.731 

.615 

2.2 

.760 

.579 

2.6 

.790 

.530 

3.0 

.818 

.474 

4.0 

.880 

.321 

5.0 

.924 

.203 

6.0 

.956 

.108 

6.5 

.968 

.069 

7.0 

.975 

.039 

7.5 

.982 

.010 

8.0 

.991 

.01 

8.5 

.996 

.01 

9.0 

.999 

.01 
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TABLE m.- FREE-STREAM DISTRIBUTIONS OF VELOCITY AND VELOCITY DERIVATIVES 
FOR ADVERSE PRESSURE GRADIENT CASE - Concluded 
[Values based on Goldberg's data (ref. 43j] 

(b) Velocity distribution b. 


dU e 



ft/ sec 

77.8 

76.7 

75.58 

74.42 

73.22 

71.98 

70.70 

69.38 

68.0 

66.56 


10.1 

.26 

65.03 

10.8 

.27 

63.43 

11.5 

.29 

61.67 

12.2 

.31 

60.75 

12.85 

.32 

58.81 

13.5 

.34 

56.88 

14.25 

.36 

55.26 

15.0 

.38 

54.03 

15.7 

.40 

52.97 

16.5 

.42 

52.05 

17.4 

.44 

51.28 

18.25 

.46 

50.65 

19.2 

.49 

50.14 

20.0 

.51 

49.74 

20.85 

.53 

49.42 

21.7 

.55 

49.16 

22.6 

.57 

49.0 

23.45 

.59 

48.9 

24.35 

.618 

48.9 

32.3 

.820 

48.9 

40.0 

1.02 

48.9 


23.7 

23.4 

23.0 

22.7 

22.3 

21.9 

21.5 

21.1 

20.7 

20.3 

19.8 

19.3 

18.8 

18.5 

17.9 

17.3 

16.8 

16.4 

16.1 

15.8 

15.6 

15.4 

15.3 
15.16 

15.06 
14.98 
14.94 

14.9 

14.9 

14.9 

14.9 


sec - - 1 

23.8 

23.9 

24.4 

24.8 

25.3 

25.5 

26.0 

26.6 

27.2 

28.4 

28.9 

31.0 

32.9 

32.7 

31.5 

25.5 

18.9 

16.0 

13.5 

11.1 

9.0 

7.18 

5.76 

4.55 

3.6 

2.54 

1.635 

.543 

0 

0 

0 


8 

8.2 

8.4 

8.6 

8.8 

9.0 

9.2 

9.4 

9.6 
9.8 

10.0 

10.2 

10.4 

10.6 

10.8 
11.0 

11.2 

11.4 

11.6 

11.8 
12.0 

12.2 

12.4 

12.6 

12.8 

13.0 

13.2 

13.4 

13.6 

15.4 

17.1 


it/ sec 

5.55X10- 5 

5.62 

5.8 
6.0 
6.2 
6.4 

6.62 

6.9 
7.21 
7.67 
8.0 
8.8 

9.6 

9.7 

9.65 
8.1 

6.15 

5.3 

4.6 
3.85 

3.15 
2.55 

2.07 

1.65 

1.3 
.93 
.6 


1.68 x 10 -5 
1.71 
1.77 
1.83 
1.89 

1.95 
2.02 
2.10 
2.20 
2.34 
2.44 

2.68 

2.93 

2.96 

2.94 
2.47 
1.87 
1.62 
1.40 
1.17 

.96 
.78 
.63 
.50 
.40 
.28 
.18 
.06 









H(R) x 2 /6>.5 

Diffusion term 

Eq. (9b) and 1. 0 

Eq. (13) 

— 1.0 

3 times eq. (13) 

— 1.0 

3 times eq. (13) 


O Wieghardt (ref. 46) 


Laminar theory (Blasius) 
C f = .664 (R Xl )‘ 1/2 

E * 
o 

2. 5 x 10' 4 
2. 5 x 10" 4 
1.0 x 10" 8 





(a) Development of velocity profiles from laminar input at R Xl = 1 x 10 4 to turbulent at R X j = 1 x 10 6 for dp/dx = 0. 

Figure 3.- Mean velocity distributions. 
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(R) x 2 /6 >.5 


Diffusion term 


Eq. (9b) 
• 1.0 
1.0 


Eq. (13) 

Eq. (13) 

3 times eq. (13) J 


R v = 1 x 10 6 
X 1 


except as noted 


Data correlation of Clauser's (ref. 42, fig. 21) 




H(R) Diffusion term 

xp/o > . 5 

Eq. (9b) Eq. (13) 

Theory ■ / — 1.0 Eq. (13) 

1.0 3 times eq. (13) 

Limits from measurements in 

equilibrium boundary layers (refs. 32 and 34) 


■ Constant value used in theory 
of Bradshaw (ref. 18) 
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X 2 /6 F=.995 


Figure 4.- Variation across boundary layer of ratio of turbulent shear stress to twice turbulent kinetic energy for flat-plate flow at R x = l x #. E 0 # = 2.5 x iq-4. 
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Experiment 



(a) Curves shown are for Eo* = 2.5 x 10”^, but for R X j = 5 x 10^, the use of Eg = 1 x 10 ® gives same curves. 

Figure 5.- The variation in /T = (u 1 ' 2 + u 2 ' 2 + across the flat-plate boundary layer as computed for various values of R X1 

with different assumptions for diffusion, and as measured by Klebanoff (ref. 32). 















Slope = 4/5 



O Wieghardt data (ref. 46) 
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Figure 9.- Comparison of calculated skin friction with measured values of Goldberg (ref. 43) for different velocity distributions, values of C, and f/S functions. 

Diffusion is 3 times Glushko diffusion and H(R) = 1 in outer part of boundary layer. 
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(ref. 43). 


2 F=. 995 

(a) xj = 12 inches (0.30 meter) and input profile used to start solution at Xj = 4 inches (0.10 meter). 

Figure 10.- Calculated distribution of \JE compared with distribution obtained from measured longitudinal turbulence intensity of Goldberg (ref. 43). 

Diffusion is 3 times Glushko diffusion. 
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(b) xj = 20 inches (0.51 meter). 
Figure 13.- Continued. 




(c) Xj = 36 inches (0.91 meter). 
Figure 13.- Concluded. 
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(a) Adverse pressure gradient region. 

Figure 16.- Comparison of the ratio of turbulent shear to turbulent kinetic energy from the theory and from the data of references 43 and 36. 
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(b) Near end of adverse pressure gradient. 

Figure 16.- Cont 
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Same as fig. 17(a) 
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(c) xj = 20 inches (0.51 meter) 


Fiaure 17.- Continued. 
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(d) Xj = 36 inches (0.91 meter). 
Figure 17.- Concluded. 
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